Study design:
#Sm_Pr$primary <- ordered(Sm_Pr$primary, levels = c("PBS", "S.m"))
table(Sm_Pr$primary, Sm_Pr$secondary_dose)
##
## 0 0.01 0.1 1 2 3 5
## PBS 115 118 206 283 263 80 165
## S.m 114 119 196 275 238 81 159
#Sm_bacterial_load$primary <- ordered(Sm_bacterial_load$primary, levels = c("PBS", "S.m"))
#table(Sm_bacterial_load$primary, Sm_bacterial_load$secondary_dose)
#THE GRAPH
#png(file="Sm_Pr.png",width=1200,height=900,res=130)
par(mar=c(4,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(Sm_Pr$day, Sm_Pr$status) ~ Sm_Pr$primary + Sm_Pr$secondary_dose), col=c("gray", "#D7EEF7", "#B5E8FC", "#09B6FC", "#055df5", "#003694", "black"), lty=c(3,3,3,3, 3, 3, 3, 1,1,1,1, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-PBS","PBS-P.rett 0.01", "PBS-P.rett 0.1", "PBS-P.rett 1.0", "PBS-P.rett 2.0","PBS-P.rett 3.0","PBS-P.rett 5.0","S.m-PBS", "S.m-P.rett 0.01", "S.m-P.rett 0.1", "S.m-P.rett 1.0", "S.m-P.rett 2.0", "S.m-P.rett 3.0", "S.m-P.rett 5.0"),col=c("gray", "#D7EEF7", "#B5E8FC", "#09B6FC", "#055df5", "#003694", "black"),lty=c(3, 3, 3, 3, 1, 1, 1, 1),lwd=5, cex=1)
Sm_Pr_PBS<-subset(Sm_Pr,secondary_dose=="0")
#table(Sm_Pr_PBS$primary)
#THE GRAPH
#png(file="AFigure1A.png",width=1200,height=900,res=130)
par(mar=c(4,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Sm_Pr_PBS), col=c("gray"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("bottomleft",inset=0, bty="n" , c("double injection control (n=115)", "chronic infection only (n=114)"),col=c("gray"),lty=c(3, 1),lwd=5, cex=1)
#dev.off()
Sm_Pr_0.01<-subset(Sm_Pr,secondary_dose=="0.01")
#table(Sm_Pr_0.01$primary)
#THE GRAPH
#png(file="Figure1A.png",width=1200,height=900,res=130)
par(mar=c(4,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Sm_Pr_0.01), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=118)", "chronic & secondary infection (n=119)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.7)
#dev.off()
Sm_Pr_0.1<-subset(Sm_Pr,secondary_dose=="0.1")
#table(Sm_Pr_0.1$primary)
#THE GRAPH
#png(file="Figure1B.png",width=1200,height=900,res=130)
par(mar=c(4,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Sm_Pr_0.1), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=206)", "chronic & secondary infection (n=196)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.7)
#dev.off()
Sm_Pr_1.0<-subset(Sm_Pr,secondary_dose=="1")
#table(Sm_Pr_1.0$primary)
#THE GRAPH
#png(file="Figure1C.png",width=1200,height=900,res=130)
par(mar=c(4,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Sm_Pr_1.0), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=283)", "chronic & secondary infection (n=275)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.7)
#dev.off()
Sm_Pr_2.0<-subset(Sm_Pr,secondary_dose=="2")
#table(Sm_Pr_2.0$primary)
#THE GRAPH
#png(file="Figure1D.png",width=1200,height=900,res=130)
par(mar=c(4,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Sm_Pr_2.0), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=263)", "chronic & secondary infection (n=238)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.7)
#dev.off()
Sm_Pr_3.0<-subset(Sm_Pr,secondary_dose=="3")
#table(Sm_Pr_3.0$primary)
#THE GRAPH
#png(file="Figure1E.png",width=1200,height=900,res=130)
par(mar=c(4,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Sm_Pr_3.0), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=80)", "chronic & secondary infection (n=81)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.7)
#dev.off()
Sm_Pr_5.0<-subset(Sm_Pr,secondary_dose=="5")
#table(Sm_Pr_5.0$primary)
#THE GRAPH
#png(file="Figure1F.png",width=1200,height=900,res=130)
par(mar=c(4,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Sm_Pr_5.0), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=165)", "chronic & secondary infection (n=159)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.7)
#dev.off()
a <- ggdraw() +
draw_image("Figure1A.png")
b <- ggdraw() +
draw_image("Figure1B.png")
c <- ggdraw() +
draw_image("Figure1C.png")
d <- ggdraw() +
draw_image("Figure1D.png")
e <- ggdraw() +
draw_image("Figure1E.png")
f <- ggdraw() +
draw_image("Figure1F.png")
multiplot <- a + b + c + d + e + f
z <- multiplot +
plot_annotation(tag_levels = 'A')
ggsave("Fig.1.ChronicSm_Pr.pdf", z, device = "pdf")
z
Sm_Pr_block1 <- subset(Sm_Pr, date == "8/2/2018")
table(Sm_Pr_block1$primary, Sm_Pr_block1$secondary_dose)
##
## 0.01 0.1 1
## PBS 17 18 18
## S.m 15 15 15
Sm_Pr_0.01_1<-subset(Sm_Pr_block1, secondary_dose=="0.01")
#table(Sm_Pr_0.01_1$primary)
Sm_Pr_0.1_1<-subset(Sm_Pr_block1, secondary_dose=="0.1")
#table(Sm_Pr_0.1_1$primary)
Sm_Pr_1.0_1<-subset(Sm_Pr_block1, secondary_dose=="1")
#table(Sm_Pr_1.0_1$primary)
#Sm_Pr_2.0_1<-subset(Sm_Pr_block1, secondary_dose=="2")
#table(Sm_Pr_2.0_1$primary)
#Sm_Pr_3.0_1<-subset(Sm_Pr_block1, secondary_dose=="3")
#table(Sm_Pr_2.0_1$primary)
#Sm_Pr_5.0_1<-subset(Sm_Pr_block1, secondary_dose=="5")
#table(Sm_Pr_5.0_1$primary)
surv_Sm_Pr_0.01_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.01_1)
surv_Sm_Pr_0.01_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.01_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 17 8 5.3 1.37 3.19
## primary=S.m 15 3 5.7 1.28 3.19
##
## Chisq= 3.2 on 1 degrees of freedom, p= 0.07
surv_Sm_Pr_0.1_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_1)
surv_Sm_Pr_0.1_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.1_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 18 18 15.5 0.395 2.33
## primary=S.m 15 12 14.5 0.423 2.33
##
## Chisq= 2.3 on 1 degrees of freedom, p= 0.1
surv_Sm_Pr_1.0_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_1)
surv_Sm_Pr_1.0_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 18 18 15.8 0.301 5.3
## primary=S.m 15 13 15.2 0.314 5.3
##
## Chisq= 5.3 on 1 degrees of freedom, p= 0.02
#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block1$day, Sm_Pr_block1$status) ~ Sm_Pr_block1$primary + Sm_Pr_block1$secondary_dose), col=c("#B5E8FC", "#09B6FC", "#055df5"), lty=c(3, 3, 3, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.01", "PBS-P.rett 0.1", "PBS-P.rett 1.0", "S.m-P.rett 0.01", "S.m-P.rett 0.1", "S.m-P.rett 1.0"),col=c("#B5E8FC", "#09B6FC", "#055df5"),lty=c(3, 3, 3, 1, 1, 1),lwd=5, cex=1)
Sm_Pr_block2 <- subset(Sm_Pr, date == "8/6/2018")
table(Sm_Pr_block2$primary, Sm_Pr_block2$secondary_dose)
##
## 0.01 0.1 1
## PBS 19 27 12
## S.m 24 23 15
Sm_Pr_0.01_2<-subset(Sm_Pr_block2, secondary_dose=="0.01")
#table(Sm_Pr_0.01_2$primary)
Sm_Pr_0.1_2<-subset(Sm_Pr_block2, secondary_dose=="0.1")
#table(Sm_Pr_0.1_2$primary)
Sm_Pr_1.0_2<-subset(Sm_Pr_block2, secondary_dose=="1")
#table(Sm_Pr_1.0_2$primary)
#Sm_Pr_2.0_1<-subset(Sm_Pr_block1, secondary_dose=="2")
#table(Sm_Pr_2.0_1$primary)
surv_Sm_Pr_0.01_2<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.01_2)
surv_Sm_Pr_0.01_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.01_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 19 13 8.24 2.74 5.57
## primary=S.m 24 9 13.76 1.64 5.57
##
## Chisq= 5.6 on 1 degrees of freedom, p= 0.02
surv_Sm_Pr_0.1_2<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_2)
surv_Sm_Pr_0.1_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.1_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 27 20 13 3.72 10.5
## primary=S.m 23 8 15 3.24 10.5
##
## Chisq= 10.5 on 1 degrees of freedom, p= 0.001
surv_Sm_Pr_1.0_2<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_2)
surv_Sm_Pr_1.0_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 12 12 8.89 1.089 4.91
## primary=S.m 15 10 13.11 0.738 4.91
##
## Chisq= 4.9 on 1 degrees of freedom, p= 0.03
#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block2$day, Sm_Pr_block2$status) ~ Sm_Pr_block2$primary + Sm_Pr_block2$secondary_dose), col=c("#B5E8FC", "#09B6FC", "#055df5"), lty=c(3, 3, 3, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.01", "PBS-P.rett 0.1", "PBS-P.rett 1.0", "S.m-P.rett 0.01", "S.m-P.rett 0.1", "S.m-P.rett 1.0"),col=c("#B5E8FC", "#09B6FC", "#055df5"),lty=c(3, 3, 3, 1, 1, 1),lwd=5, cex=1)
Sm_Pr_block3 <- subset(Sm_Pr, date == "6/14/2019")
table(Sm_Pr_block3$primary, Sm_Pr_block3$secondary_dose)
##
## 0.01 0.1 1 2
## PBS 15 14 17 20
## S.m 15 19 17 18
Sm_Pr_0.01_3<-subset(Sm_Pr_block3, secondary_dose=="0.01")
#table(Sm_Pr_0.01_3$primary)
Sm_Pr_0.1_3<-subset(Sm_Pr_block3, secondary_dose=="0.1")
#table(Sm_Pr_0.1_3$primary)
Sm_Pr_1.0_3<-subset(Sm_Pr_block3, secondary_dose=="1")
#table(Sm_Pr_1.0_3$primary)
Sm_Pr_2.0_3<-subset(Sm_Pr_block3, secondary_dose=="2")
#table(Sm_Pr_2.0_3$primary)
surv_Sm_Pr_0.01_3<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.01_3)
surv_Sm_Pr_0.01_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.01_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 15 9 6.24 1.22 3.41
## primary=S.m 15 4 6.76 1.13 3.41
##
## Chisq= 3.4 on 1 degrees of freedom, p= 0.06
surv_Sm_Pr_0.1_3<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_3)
surv_Sm_Pr_0.1_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.1_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 14 7 2.74 6.63 12.1
## primary=S.m 19 0 4.26 4.26 12.1
##
## Chisq= 12.1 on 1 degrees of freedom, p= 5e-04
surv_Sm_Pr_1.0_3<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_3)
surv_Sm_Pr_1.0_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 17 17 8.5 8.50 33
## primary=S.m 17 3 11.5 6.28 33
##
## Chisq= 33 on 1 degrees of freedom, p= 9e-09
surv_Sm_Pr_2.0_3<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_2.0_3)
surv_Sm_Pr_2.0_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_2.0_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 20 20 12.1 5.15 26.8
## primary=S.m 18 5 12.9 4.83 26.8
##
## Chisq= 26.8 on 1 degrees of freedom, p= 2e-07
#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block3$day, Sm_Pr_block3$status) ~ Sm_Pr_block3$primary + Sm_Pr_block3$secondary_dose), col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 3, 1, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.01", "PBS-P.rett 0.1", "PBS-P.rett 1.0", "PBS-P.rett 2.0", "S.m-P.rett 0.01", "S.m-P.rett 0.1", "S.m-P.rett 1.0", "S.m-P.rett 2.0"),col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 3, 1, 1, 1, 1),lwd=5, cex=1)
Sm_Pr_block4 <- subset(Sm_Pr, date == "6/17/2019")
table(Sm_Pr_block4$primary, Sm_Pr_block4$secondary_dose)
##
## 0.01 0.1 1 2
## PBS 15 19 17 17
## S.m 18 13 16 11
Sm_Pr_0.01_4<-subset(Sm_Pr_block4, secondary_dose=="0.01")
#table(Sm_Pr_0.01_4$primary)
Sm_Pr_0.1_4<-subset(Sm_Pr_block4, secondary_dose=="0.1")
#table(Sm_Pr_0.1_4$primary)
Sm_Pr_1.0_4<-subset(Sm_Pr_block4, secondary_dose=="1")
#table(Sm_Pr_1.0_4$primary)
Sm_Pr_2.0_4<-subset(Sm_Pr_block4, secondary_dose=="2")
#table(Sm_Pr_2.0_4$primary)
surv_Sm_Pr_0.01_4<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.01_4)
surv_Sm_Pr_0.01_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.01_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 15 10 7.59 0.767 1.76
## primary=S.m 18 8 10.41 0.559 1.76
##
## Chisq= 1.8 on 1 degrees of freedom, p= 0.2
surv_Sm_Pr_0.1_4<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_4)
surv_Sm_Pr_0.1_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.1_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 19 9 7.85 0.170 0.438
## primary=S.m 13 5 6.15 0.216 0.438
##
## Chisq= 0.4 on 1 degrees of freedom, p= 0.5
surv_Sm_Pr_1.0_4<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_4)
surv_Sm_Pr_1.0_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 17 16 9.76 4.00 15.7
## primary=S.m 16 6 12.24 3.18 15.7
##
## Chisq= 15.7 on 1 degrees of freedom, p= 7e-05
surv_Sm_Pr_2.0_4<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_2.0_4)
surv_Sm_Pr_2.0_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_2.0_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 17 17 12.14 1.94 16.7
## primary=S.m 11 5 9.86 2.39 16.7
##
## Chisq= 16.7 on 1 degrees of freedom, p= 4e-05
#surv_Sm_Pr_3.0_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_3.0_1)
#surv_Sm_Pr_3.0_1
#surv_Sm_Pr_5.0_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_5.0_1)
#surv_Sm_Pr_5.0_1
#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block4$day, Sm_Pr_block4$status) ~ Sm_Pr_block4$primary + Sm_Pr_block4$secondary_dose), col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 3, 1, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.01", "PBS-P.rett 0.1", "PBS-P.rett 1.0", "PBS-P.rett 2.0", "S.m-P.rett 0.01", "S.m-P.rett 0.1", "S.m-P.rett 1.0", "S.m-P.rett 2.0"),col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 3, 1, 1, 1, 1),lwd=5, cex=1)
Sm_Pr_block5 <- subset(Sm_Pr, date == "6/19/2019")
table(Sm_Pr_block5$primary, Sm_Pr_block5$secondary_dose)
##
## 0.01 0.1 1 2
## PBS 15 5 9 19
## S.m 11 17 14 15
Sm_Pr_0.01_5<-subset(Sm_Pr_block5, secondary_dose=="0.01")
#table(Sm_Pr_0.01_5$primary)
Sm_Pr_0.1_5<-subset(Sm_Pr_block5, secondary_dose=="0.1")
#table(Sm_Pr_0.1_5$primary)
Sm_Pr_1.0_5<-subset(Sm_Pr_block5, secondary_dose=="1")
#table(Sm_Pr_1.0_5$primary)
Sm_Pr_2.0_5<-subset(Sm_Pr_block5, secondary_dose=="2")
#table(Sm_Pr_2.0_5$primary)
surv_Sm_Pr_0.01_5<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.01_5)
surv_Sm_Pr_0.01_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.01_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 15 10 7.83 0.604 1.61
## primary=S.m 11 5 7.17 0.659 1.61
##
## Chisq= 1.6 on 1 degrees of freedom, p= 0.2
surv_Sm_Pr_0.1_5<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_5)
surv_Sm_Pr_0.1_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.1_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 5 5 1.21 11.90 17.1
## primary=S.m 17 5 8.79 1.64 17.1
##
## Chisq= 17.1 on 1 degrees of freedom, p= 4e-05
surv_Sm_Pr_1.0_5<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_5)
surv_Sm_Pr_1.0_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 9 9 5.29 2.61 8.8
## primary=S.m 14 7 10.71 1.29 8.8
##
## Chisq= 8.8 on 1 degrees of freedom, p= 0.003
surv_Sm_Pr_2.0_5<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_2.0_5)
surv_Sm_Pr_2.0_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_2.0_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 19 18 11.8 3.27 15
## primary=S.m 15 9 15.2 2.54 15
##
## Chisq= 15 on 1 degrees of freedom, p= 1e-04
#surv_Sm_Pr_3.0_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_3.0_1)
#surv_Sm_Pr_3.0_1
#surv_Sm_Pr_5.0_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_5.0_1)
#surv_Sm_Pr_5.0_1
#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block5$day, Sm_Pr_block5$status) ~ Sm_Pr_block5$primary + Sm_Pr_block5$secondary_dose), col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 3, 1, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.01", "PBS-P.rett 0.1", "PBS-P.rett 1.0", "PBS-P.rett 2.0", "S.m-P.rett 0.01", "S.m-P.rett 0.1", "S.m-P.rett 1.0", "S.m-P.rett 2.0"),col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 3, 1, 1, 1, 1),lwd=5, cex=1)
Sm_Pr_block6 <- subset(Sm_Pr, date == "6/25/2019")
table(Sm_Pr_block6$primary, Sm_Pr_block6$secondary_dose)
##
## 0.01 0.1 1 2
## PBS 17 19 19 20
## S.m 20 15 17 16
Sm_Pr_0.01_6<-subset(Sm_Pr_block6, secondary_dose=="0.01")
#table(Sm_Pr_0.01_6$primary)
Sm_Pr_0.1_6<-subset(Sm_Pr_block6, secondary_dose=="0.1")
#table(Sm_Pr_0.1_6$primary)
Sm_Pr_1.0_6<-subset(Sm_Pr_block6, secondary_dose=="1")
#table(Sm_Pr_1.0_6$primary)
Sm_Pr_2.0_6<-subset(Sm_Pr_block6, secondary_dose=="2")
#table(Sm_Pr_2.0_6$primary)
surv_Sm_Pr_0.01_6<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.01_6)
surv_Sm_Pr_0.01_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.01_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 17 12 7.11 3.36 7.48
## primary=S.m 20 5 9.89 2.42 7.48
##
## Chisq= 7.5 on 1 degrees of freedom, p= 0.006
surv_Sm_Pr_0.1_6<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_6)
surv_Sm_Pr_0.1_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.1_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 19 15 8.52 4.92 13.2
## primary=S.m 15 2 8.48 4.95 13.2
##
## Chisq= 13.2 on 1 degrees of freedom, p= 3e-04
surv_Sm_Pr_1.0_6<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_6)
surv_Sm_Pr_1.0_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 19 17 9.99 4.92 18
## primary=S.m 17 5 12.01 4.09 18
##
## Chisq= 18 on 1 degrees of freedom, p= 2e-05
surv_Sm_Pr_2.0_6<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_2.0_6)
surv_Sm_Pr_2.0_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_2.0_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 20 20 11.9 5.57 27.8
## primary=S.m 16 6 14.1 4.68 27.8
##
## Chisq= 27.8 on 1 degrees of freedom, p= 1e-07
#surv_Sm_Pr_3.0_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_3.0_1)
#surv_Sm_Pr_3.0_1
#surv_Sm_Pr_5.0_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_5.0_1)
#surv_Sm_Pr_5.0_1
#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block6$day, Sm_Pr_block6$status) ~ Sm_Pr_block6$primary + Sm_Pr_block6$secondary_dose), col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 3, 1, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.01", "PBS-P.rett 0.1", "PBS-P.rett 1.0", "PBS-P.rett 2.0", "S.m-P.rett 0.01", "S.m-P.rett 0.1", "S.m-P.rett 1.0", "S.m-P.rett 2.0"),col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 3, 1, 1, 1, 1),lwd=5, cex=1)
Sm_Pr_block7 <- subset(Sm_Pr, date == "6/27/2019")
table(Sm_Pr_block7$primary, Sm_Pr_block7$secondary_dose)
##
## 0.01 0.1 1 2
## PBS 20 17 20 19
## S.m 16 14 14 14
Sm_Pr_0.01_7<-subset(Sm_Pr_block7, secondary_dose=="0.01")
#table(Sm_Pr_0.01_7$primary)
Sm_Pr_0.1_7<-subset(Sm_Pr_block7, secondary_dose=="0.1")
#table(Sm_Pr_0.1_7$primary)
Sm_Pr_1.0_7<-subset(Sm_Pr_block7, secondary_dose=="1")
#table(Sm_Pr_1.0_7$primary)
Sm_Pr_2.0_7<-subset(Sm_Pr_block7, secondary_dose=="2")
#table(Sm_Pr_2.0_7$primary)
surv_Sm_Pr_0.01_7<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.01_7)
surv_Sm_Pr_0.01_7
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.01_7)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 20 11 10.11 0.0779 0.206
## primary=S.m 16 9 9.89 0.0797 0.206
##
## Chisq= 0.2 on 1 degrees of freedom, p= 0.7
surv_Sm_Pr_0.1_7<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_7)
surv_Sm_Pr_0.1_7
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.1_7)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 17 14 9.38 2.27 6.94
## primary=S.m 14 6 10.62 2.01 6.94
##
## Chisq= 6.9 on 1 degrees of freedom, p= 0.008
surv_Sm_Pr_1.0_7<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_7)
surv_Sm_Pr_1.0_7
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_7)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 20 20 13.5 3.09 22.5
## primary=S.m 14 7 13.5 3.11 22.5
##
## Chisq= 22.5 on 1 degrees of freedom, p= 2e-06
surv_Sm_Pr_2.0_7<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_2.0_7)
surv_Sm_Pr_2.0_7
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_2.0_7)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 19 18 11.2 4.18 17.4
## primary=S.m 14 9 15.8 2.95 17.4
##
## Chisq= 17.4 on 1 degrees of freedom, p= 3e-05
#surv_Sm_Pr_3.0_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_3.0_1)
#surv_Sm_Pr_3.0_1
#surv_Sm_Pr_5.0_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_5.0_1)
#surv_Sm_Pr_5.0_1
#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block7$day, Sm_Pr_block7$status) ~ Sm_Pr_block7$primary + Sm_Pr_block7$secondary_dose), col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 3, 1, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.01", "PBS-P.rett 0.1", "PBS-P.rett 1.0", "PBS-P.rett 2.0", "S.m-P.rett 0.01", "S.m-P.rett 0.1", "S.m-P.rett 1.0", "S.m-P.rett 2.0"),col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 3, 1, 1, 1, 1),lwd=5, cex=1)
Sm_Pr_block8 <- subset(Sm_Pr, date == "1/30/2020")
table(Sm_Pr_block8$primary, Sm_Pr_block8$secondary_dose)
##
## 0 0.1 1 2 5
## PBS 18 29 27 26 28
## S.m 20 23 24 25 25
Sm_Pr_0.1_8<-subset(Sm_Pr_block8, secondary_dose=="0.1")
#table(Sm_Pr_0.1_8$primary)
Sm_Pr_1.0_8<-subset(Sm_Pr_block8, secondary_dose=="1")
#table(Sm_Pr_1.0_8$primary)
Sm_Pr_2.0_8<-subset(Sm_Pr_block8, secondary_dose=="2")
#table(Sm_Pr_2.0_8$primary)
#Sm_Pr_3.0_8<-subset(Sm_Pr_block8, secondary_dose=="3")
#table(Sm_Pr_3.0_8$primary)
Sm_Pr_5.0_8<-subset(Sm_Pr_block8, secondary_dose=="5")
#table(Sm_Pr_5.0_8$primary)
surv_Sm_Pr_0.1_8<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_8)
surv_Sm_Pr_0.1_8
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.1_8)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 29 29 15.7 11.33 26.8
## primary=S.m 23 15 28.3 6.27 26.8
##
## Chisq= 26.8 on 1 degrees of freedom, p= 2e-07
surv_Sm_Pr_1.0_8<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_8)
surv_Sm_Pr_1.0_8
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_8)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 27 27 16.9 5.97 33.4
## primary=S.m 24 19 29.1 3.48 33.4
##
## Chisq= 33.4 on 1 degrees of freedom, p= 8e-09
surv_Sm_Pr_2.0_8<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_2.0_8)
surv_Sm_Pr_2.0_8
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_2.0_8)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 26 26 15.3 7.49 36.4
## primary=S.m 25 20 30.7 3.73 36.4
##
## Chisq= 36.4 on 1 degrees of freedom, p= 2e-09
#surv_Sm_Pr_3.0_8<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_3.0_8)
#surv_Sm_Pr_3.0_8
surv_Sm_Pr_5.0_8<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_5.0_8)
surv_Sm_Pr_5.0_8
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_5.0_8)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 28 28 15.3 10.49 48.2
## primary=S.m 25 16 28.7 5.61 48.2
##
## Chisq= 48.2 on 1 degrees of freedom, p= 4e-12
#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block8$day, Sm_Pr_block8$status) ~ Sm_Pr_block8$primary + Sm_Pr_block8$secondary_dose), col=c("gray", "#B5E8FC", "#09B6FC", "#055df5", "black"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1), lwd=5, yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-PBS", "PBS-P.rett 0.1", "PBS-P.rett 1.0", "PBS-P.rett 2.0", "PBS-P.rett 5.0", "S.m-PBS", "S.m-P.rett 0.1", "S.m-P.rett 1.0", "S.m-P.rett 2.0", "S.m-P.rett 5.0"),col=c("gray", "#B5E8FC", "#09B6FC", "#055df5", "black"),lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1),lwd=5, cex=1)
Sm_Pr_block9 <- subset(Sm_Pr, date == "2/13/2020")
table(Sm_Pr_block9$primary, Sm_Pr_block9$secondary_dose)
##
## 0 0.1 1 2 5
## PBS 18 28 31 30 30
## S.m 18 28 28 28 28
Sm_Pr_0.1_9<-subset(Sm_Pr_block9, secondary_dose=="0.1")
#table(Sm_Pr_0.1_9$primary)
Sm_Pr_1.0_9<-subset(Sm_Pr_block9, secondary_dose=="1")
#table(Sm_Pr_1.0_9$primary)
Sm_Pr_2.0_9<-subset(Sm_Pr_block9, secondary_dose=="2")
#table(Sm_Pr_2.0_9$primary)
#Sm_Pr_3.0_9<-subset(Sm_Pr_block9, secondary_dose=="3")
#table(Sm_Pr_3.0_9$primary)
Sm_Pr_5.0_9<-subset(Sm_Pr_block9, secondary_dose=="5")
#table(Sm_Pr_5.0_9$primary)
surv_Sm_Pr_0.1_9<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_9)
surv_Sm_Pr_0.1_9
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.1_9)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 28 26 13.9 10.64 21.8
## primary=S.m 28 10 22.1 6.66 21.8
##
## Chisq= 21.8 on 1 degrees of freedom, p= 3e-06
surv_Sm_Pr_1.0_9<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_9)
surv_Sm_Pr_1.0_9
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_9)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 31 31 18.9 7.72 41
## primary=S.m 28 15 27.1 5.39 41
##
## Chisq= 41 on 1 degrees of freedom, p= 2e-10
surv_Sm_Pr_2.0_9<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_2.0_9)
surv_Sm_Pr_2.0_9
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_2.0_9)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30 30 17.1 9.80 46.3
## primary=S.m 28 8 20.9 7.99 46.3
##
## Chisq= 46.3 on 1 degrees of freedom, p= 1e-11
#surv_Sm_Pr_3.0_9<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_3.0_9)
#surv_Sm_Pr_3.0_9
surv_Sm_Pr_5.0_9<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_5.0_9)
surv_Sm_Pr_5.0_9
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_5.0_9)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30 30 18.2 7.69 35.3
## primary=S.m 28 21 32.8 4.26 35.3
##
## Chisq= 35.3 on 1 degrees of freedom, p= 3e-09
#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block9$day, Sm_Pr_block9$status) ~ Sm_Pr_block9$primary + Sm_Pr_block9$secondary_dose), col=c("gray", "#B5E8FC", "#09B6FC", "#055df5", "black"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1), lwd=5, yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-PBS", "PBS-P.rett 0.1", "PBS-P.rett 1.0", "PBS-P.rett 2.0", "PBS-P.rett 5.0", "S.m-PBS", "S.m-P.rett 0.1", "S.m-P.rett 1.0", "S.m-P.rett 2.0", "S.m-P.rett 5.0"),col=c("gray", "#B5E8FC", "#09B6FC", "#055df5", "black"),lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1),lwd=5, cex=1)
Sm_Pr_block10 <- subset(Sm_Pr, date == "9/27/2020")
table(Sm_Pr_block10$primary, Sm_Pr_block10$secondary_dose)
##
## 0 1 2 3 5
## PBS 10 27 26 20 20
## S.m 10 28 27 28 27
#Sm_Pr_0.1_10<-subset(Sm_Pr_block10, secondary_dose=="0.1")
#table(Sm_Pr_0.1_10$primary)
Sm_Pr_1.0_10<-subset(Sm_Pr_block10, secondary_dose=="1")
#table(Sm_Pr_1.0_10$primary)
Sm_Pr_2.0_10<-subset(Sm_Pr_block10, secondary_dose=="2")
#table(Sm_Pr_2.0_10$primary)
Sm_Pr_3.0_10<-subset(Sm_Pr_block10, secondary_dose=="3")
#table(Sm_Pr_3.0_10$primary)
Sm_Pr_5.0_10<-subset(Sm_Pr_block10, secondary_dose=="5")
#table(Sm_Pr_5.0_10$primary)
#surv_Sm_Pr_0.1_10<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_10)
#surv_Sm_Pr_0.1_10
surv_Sm_Pr_1.0_10<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_10)
surv_Sm_Pr_1.0_10
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_10)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 27 27 12.7 16.16 53.6
## primary=S.m 28 7 21.3 9.62 53.6
##
## Chisq= 53.6 on 1 degrees of freedom, p= 2e-13
surv_Sm_Pr_2.0_10<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_2.0_10)
surv_Sm_Pr_2.0_10
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_2.0_10)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 26 26 14.7 8.65 38.4
## primary=S.m 27 9 20.3 6.28 38.4
##
## Chisq= 38.4 on 1 degrees of freedom, p= 6e-10
surv_Sm_Pr_3.0_10<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_3.0_10)
surv_Sm_Pr_3.0_10
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_3.0_10)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 20 20 10.4 8.82 30.9
## primary=S.m 28 10 19.6 4.69 30.9
##
## Chisq= 30.9 on 1 degrees of freedom, p= 3e-08
surv_Sm_Pr_5.0_10<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_5.0_10)
surv_Sm_Pr_5.0_10
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_5.0_10)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 20 20 12.8 4.10 19.3
## primary=S.m 27 16 23.2 2.25 19.3
##
## Chisq= 19.3 on 1 degrees of freedom, p= 1e-05
#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block10$day, Sm_Pr_block10$status) ~ Sm_Pr_block10$primary + Sm_Pr_block10$secondary_dose), col=c("gray", "#09B6FC", "#055df5", "#003694", "black"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1), lwd=5, yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-PBS", "PBS-P.rett 1.0", "PBS-P.rett 2.0", "PBS-P.rett 3.0", "PBS-P.rett 5.0", "S.m-PBS", "S.m-P.rett 1.0", "S.m-P.rett 2.0", "S.m-P.rett 3.0", "S.m-P.rett 5.0"),col=c("gray", "#09B6FC", "#055df5", "#003694", "black"),lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1),lwd=5, cex=1)
Sm_Pr_block11 <- subset(Sm_Pr, date == "10/6/2020")
table(Sm_Pr_block11$primary, Sm_Pr_block11$secondary_dose)
##
## 0 1 2 3 5
## PBS 20 28 28 30 29
## S.m 16 28 27 26 27
#Sm_Pr_0.1_11<-subset(Sm_Pr_block11, secondary_dose=="0.1")
#table(Sm_Pr_0.1_11$primary)
Sm_Pr_1.0_11<-subset(Sm_Pr_block11, secondary_dose=="1")
#table(Sm_Pr_1.0_11$primary)
Sm_Pr_2.0_11<-subset(Sm_Pr_block11, secondary_dose=="2")
#table(Sm_Pr_2.0_11$primary)
Sm_Pr_3.0_11<-subset(Sm_Pr_block11, secondary_dose=="3")
#table(Sm_Pr_3.0_11$primary)
Sm_Pr_5.0_11<-subset(Sm_Pr_block11, secondary_dose=="5")
#table(Sm_Pr_5.0_11$primary)
#surv_Sm_Pr_0.1_11<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_11)
#surv_Sm_Pr_0.1_11
surv_Sm_Pr_1.0_11<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_11)
surv_Sm_Pr_1.0_11
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_11)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 28 24 14.2 6.85 19
## primary=S.m 28 12 21.8 4.44 19
##
## Chisq= 19 on 1 degrees of freedom, p= 1e-05
surv_Sm_Pr_2.0_11<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_2.0_11)
surv_Sm_Pr_2.0_11
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_2.0_11)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 28 28 14.4 12.93 51.5
## primary=S.m 27 9 22.6 8.21 51.5
##
## Chisq= 51.5 on 1 degrees of freedom, p= 7e-13
surv_Sm_Pr_3.0_11<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_3.0_11)
surv_Sm_Pr_3.0_11
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_3.0_11)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30 30 16.1 12.07 55
## primary=S.m 26 14 27.9 6.95 55
##
## Chisq= 55 on 1 degrees of freedom, p= 1e-13
surv_Sm_Pr_5.0_11<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_5.0_11)
surv_Sm_Pr_5.0_11
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_5.0_11)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 29 29 16.1 10.4 47.6
## primary=S.m 27 11 23.9 7.0 47.6
##
## Chisq= 47.6 on 1 degrees of freedom, p= 5e-12
#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block11$day, Sm_Pr_block11$status) ~ Sm_Pr_block11$primary + Sm_Pr_block11$secondary_dose), col=c("gray", "#09B6FC", "#055df5", "#003694", "black"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1), lwd=5, yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-PBS", "PBS-P.rett 1.0", "PBS-P.rett 2.0", "PBS-P.rett 3.0", "PBS-P.rett 5.0", "S.m-PBS", "S.m-P.rett 1.0", "S.m-P.rett 2.0", "S.m-P.rett 3.0", "S.m-P.rett 5.0"),col=c("gray", "#09B6FC", "#055df5", "#003694", "black"),lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1),lwd=5, cex=1)
Sm_Pr_block12 <- subset(Sm_Pr, date == "11/1/2020")
table(Sm_Pr_block12$primary, Sm_Pr_block12$secondary_dose)
##
## 0 1 2 3 5
## PBS 20 30 30 30 30
## S.m 20 30 30 27 25
#Sm_Pr_0.1_12<-subset(Sm_Pr_block12, secondary_dose=="0.1")
#table(Sm_Pr_0.1_12$primary)
Sm_Pr_1.0_12<-subset(Sm_Pr_block12, secondary_dose=="1")
#table(Sm_Pr_1.0_12$primary)
Sm_Pr_2.0_12<-subset(Sm_Pr_block12, secondary_dose=="2")
#table(Sm_Pr_2.0_12$primary)
Sm_Pr_3.0_12<-subset(Sm_Pr_block12, secondary_dose=="3")
#table(Sm_Pr_3.0_12$primary)
Sm_Pr_5.0_12<-subset(Sm_Pr_block12, secondary_dose=="5")
#table(Sm_Pr_5.0_12$primary)
#surv_Sm_Pr_0.1_12<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_12)
#surv_Sm_Pr_0.1_12
surv_Sm_Pr_1.0_12<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_12)
surv_Sm_Pr_1.0_12
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_12)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30 28 13.6 15.17 33.2
## primary=S.m 30 13 27.4 7.55 33.2
##
## Chisq= 33.2 on 1 degrees of freedom, p= 8e-09
surv_Sm_Pr_2.0_12<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_2.0_12)
surv_Sm_Pr_2.0_12
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_2.0_12)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30 30 13.9 18.5 58.6
## primary=S.m 30 10 26.1 9.9 58.6
##
## Chisq= 58.6 on 1 degrees of freedom, p= 2e-14
surv_Sm_Pr_3.0_12<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_3.0_12)
surv_Sm_Pr_3.0_12
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_3.0_12)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30 30 15.6 13.39 50
## primary=S.m 27 21 35.4 5.88 50
##
## Chisq= 50 on 1 degrees of freedom, p= 2e-12
surv_Sm_Pr_5.0_12<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_5.0_12)
surv_Sm_Pr_5.0_12
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_5.0_12)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30 30 20.2 4.78 31.5
## primary=S.m 25 17 26.8 3.59 31.5
##
## Chisq= 31.5 on 1 degrees of freedom, p= 2e-08
#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block12$day, Sm_Pr_block12$status) ~ Sm_Pr_block12$primary + Sm_Pr_block12$secondary_dose), col=c("gray", "#09B6FC", "#055df5", "#003694", "black"), lty=c(3, 3, 3, 3, 3, 3, 1, 1, 1, 1, 1, 1), lwd=5, yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-PBS", "PBS-P.rett 1.0", "PBS-P.rett 2.0", "PBS-P.rett 3.0", "PBS-P.rett 5.0", "S.m-PBS", "S.m-P.rett 1.0", "S.m-P.rett 2.0", "S.m-P.rett 3.0", "S.m-P.rett 5.0"),col=c("gray", "#09B6FC", "#055df5", "#003694", "black"),lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1),lwd=5, cex=1)
Sm_Pr_block13 <- subset(Sm_Pr, date == "6/16/2021")
table(Sm_Pr_block13$primary, Sm_Pr_block13$secondary_dose)
##
## 0 0.1 1 2 5
## PBS 29 30 28 28 28
## S.m 30 29 29 27 27
Sm_Pr_0.1_13<-subset(Sm_Pr_block13, secondary_dose=="0.1")
#table(Sm_Pr_0.1_13$primary)
Sm_Pr_1.0_13<-subset(Sm_Pr_block13, secondary_dose=="1")
#table(Sm_Pr_1.0_13$primary)
Sm_Pr_2.0_13<-subset(Sm_Pr_block13, secondary_dose=="2")
#table(Sm_Pr_2.0_13$primary)
#Sm_Pr_3.0_13<-subset(Sm_Pr_block13, secondary_dose=="3")
#table(Sm_Pr_2.0_13$primary)
Sm_Pr_5.0_13<-subset(Sm_Pr_block13, secondary_dose=="5")
#table(Sm_Pr_2.0_13$primary)
surv_Sm_Pr_0.1_13<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_13)
surv_Sm_Pr_0.1_13
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.1_13)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30 28 14.5 12.56 39.3
## primary=S.m 29 6 19.5 9.34 39.3
##
## Chisq= 39.3 on 1 degrees of freedom, p= 4e-10
surv_Sm_Pr_1.0_13<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_13)
surv_Sm_Pr_1.0_13
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_13)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 28 28 14.3 13.2 51.4
## primary=S.m 29 4 17.7 10.6 51.4
##
## Chisq= 51.4 on 1 degrees of freedom, p= 7e-13
surv_Sm_Pr_2.0_13<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_2.0_13)
surv_Sm_Pr_2.0_13
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_2.0_13)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 28 28 15.4 10.42 45.1
## primary=S.m 27 5 17.6 9.06 45.1
##
## Chisq= 45.1 on 1 degrees of freedom, p= 2e-11
#surv_Sm_Pr_3.0_13<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_3.0_13)
#surv_Sm_Pr_3.0_13
surv_Sm_Pr_5.0_13<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_5.0_13)
surv_Sm_Pr_5.0_13
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_5.0_13)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 28 28 15.4 10.42 45.1
## primary=S.m 27 4 16.6 9.61 45.1
##
## Chisq= 45.1 on 1 degrees of freedom, p= 2e-11
#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block13$day, Sm_Pr_block13$status) ~ Sm_Pr_block13$primary + Sm_Pr_block13$secondary_dose), col=c("gray", "#09B6FC", "#055df5", "#003694", "black"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1), lwd=5, yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-PBS", "PBS-P.rett 1.0", "PBS-P.rett 2.0", "PBS-P.rett 3.0", "PBS-P.rett 5.0", "S.m-PBS", "S.m-P.rett 1.0", "S.m-P.rett 2.0", "S.m-P.rett 3.0", "S.m-P.rett 5.0"),col=c("gray", "#09B6FC", "#055df5", "#003694", "black"),lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1),lwd=5, cex=1)
# Grab p-values from each individual logrank test (results shown above in each block)
pvalue_0.01_1<-pchisq(surv_Sm_Pr_0.01_1$chisq, df=1, lower.tail=F)
pvalue_0.1_1<-pchisq(surv_Sm_Pr_0.1_1$chisq, df=1, lower.tail=F)
pvalue_1.0_1<-pchisq(surv_Sm_Pr_1.0_1$chisq, df=1, lower.tail=F)
pvalue_0.01_2<-pchisq(surv_Sm_Pr_0.01_2$chisq, df=1, lower.tail=F)
pvalue_0.1_2<-pchisq(surv_Sm_Pr_0.1_2$chisq, df=1, lower.tail=F)
pvalue_1.0_2<-pchisq(surv_Sm_Pr_1.0_2$chisq, df=1, lower.tail=F)
pvalue_0.01_3<-pchisq(surv_Sm_Pr_0.01_3$chisq, df=1, lower.tail=F)
pvalue_0.1_3<-pchisq(surv_Sm_Pr_0.1_3$chisq, df=1, lower.tail=F)
pvalue_1.0_3<-pchisq(surv_Sm_Pr_1.0_3$chisq, df=1, lower.tail=F)
pvalue_2.0_3<-pchisq(surv_Sm_Pr_2.0_3$chisq, df=1, lower.tail=F)
pvalue_0.01_4<-pchisq(surv_Sm_Pr_0.01_4$chisq, df=1, lower.tail=F)
pvalue_0.1_4<-pchisq(surv_Sm_Pr_0.1_4$chisq, df=1, lower.tail=F)
pvalue_1.0_4<-pchisq(surv_Sm_Pr_1.0_4$chisq, df=1, lower.tail=F)
pvalue_2.0_4<-pchisq(surv_Sm_Pr_2.0_4$chisq, df=1, lower.tail=F)
pvalue_0.01_5<-pchisq(surv_Sm_Pr_0.01_5$chisq, df=1, lower.tail=F)
pvalue_0.1_5<-pchisq(surv_Sm_Pr_0.1_5$chisq, df=1, lower.tail=F)
pvalue_1.0_5<-pchisq(surv_Sm_Pr_1.0_5$chisq, df=1, lower.tail=F)
pvalue_2.0_5<-pchisq(surv_Sm_Pr_2.0_5$chisq, df=1, lower.tail=F)
pvalue_0.01_6<-pchisq(surv_Sm_Pr_0.01_6$chisq, df=1, lower.tail=F)
pvalue_0.1_6<-pchisq(surv_Sm_Pr_0.1_6$chisq, df=1, lower.tail=F)
pvalue_1.0_6<-pchisq(surv_Sm_Pr_1.0_6$chisq, df=1, lower.tail=F)
pvalue_2.0_6<-pchisq(surv_Sm_Pr_2.0_6$chisq, df=1, lower.tail=F)
pvalue_0.01_7<-pchisq(surv_Sm_Pr_0.01_7$chisq, df=1, lower.tail=F)
pvalue_0.1_7<-pchisq(surv_Sm_Pr_0.1_7$chisq, df=1, lower.tail=F)
pvalue_1.0_7<-pchisq(surv_Sm_Pr_1.0_7$chisq, df=1, lower.tail=F)
pvalue_2.0_7<-pchisq(surv_Sm_Pr_2.0_7$chisq, df=1, lower.tail=F)
pvalue_0.1_8<-pchisq(surv_Sm_Pr_0.1_8$chisq, df=1, lower.tail=F)
pvalue_1.0_8<-pchisq(surv_Sm_Pr_1.0_8$chisq, df=1, lower.tail=F)
pvalue_2.0_8<-pchisq(surv_Sm_Pr_2.0_8$chisq, df=1, lower.tail=F)
pvalue_5.0_8<-pchisq(surv_Sm_Pr_5.0_8$chisq, df=1, lower.tail=F)
pvalue_0.1_9<-pchisq(surv_Sm_Pr_0.1_9$chisq, df=1, lower.tail=F)
pvalue_1.0_9<-pchisq(surv_Sm_Pr_1.0_9$chisq, df=1, lower.tail=F)
pvalue_2.0_9<-pchisq(surv_Sm_Pr_2.0_9$chisq, df=1, lower.tail=F)
pvalue_5.0_9<-pchisq(surv_Sm_Pr_5.0_9$chisq, df=1, lower.tail=F)
pvalue_1.0_10<-pchisq(surv_Sm_Pr_1.0_10$chisq, df=1, lower.tail=F)
pvalue_2.0_10<-pchisq(surv_Sm_Pr_2.0_10$chisq, df=1, lower.tail=F)
pvalue_3.0_10<-pchisq(surv_Sm_Pr_3.0_10$chisq, df=1, lower.tail=F)
pvalue_5.0_10<-pchisq(surv_Sm_Pr_5.0_10$chisq, df=1, lower.tail=F)
pvalue_1.0_11<-pchisq(surv_Sm_Pr_1.0_11$chisq, df=1, lower.tail=F)
pvalue_2.0_11<-pchisq(surv_Sm_Pr_2.0_11$chisq, df=1, lower.tail=F)
pvalue_3.0_11<-pchisq(surv_Sm_Pr_3.0_11$chisq, df=1, lower.tail=F)
pvalue_5.0_11<-pchisq(surv_Sm_Pr_5.0_11$chisq, df=1, lower.tail=F)
pvalue_1.0_12<-pchisq(surv_Sm_Pr_1.0_12$chisq, df=1, lower.tail=F)
pvalue_2.0_12<-pchisq(surv_Sm_Pr_2.0_12$chisq, df=1, lower.tail=F)
pvalue_3.0_12<-pchisq(surv_Sm_Pr_3.0_12$chisq, df=1, lower.tail=F)
pvalue_5.0_12<-pchisq(surv_Sm_Pr_5.0_12$chisq, df=1, lower.tail=F)
pvalue_0.1_13<-pchisq(surv_Sm_Pr_0.1_13$chisq, df=1, lower.tail=F)
pvalue_1.0_13<-pchisq(surv_Sm_Pr_1.0_13$chisq, df=1, lower.tail=F)
pvalue_2.0_13<-pchisq(surv_Sm_Pr_2.0_13$chisq, df=1, lower.tail=F)
pvalue_5.0_13<-pchisq(surv_Sm_Pr_5.0_13$chisq, df=1, lower.tail=F)
# Use all p-values from a given condition to get Fisher's combined test p-value
pvalues_0.01 <- c(pvalue_0.01_1,pvalue_0.01_2,pvalue_0.01_3,pvalue_0.01_4,pvalue_0.01_5,pvalue_0.01_6,pvalue_0.01_7)
test.stat_0.01 <- -2*sum(log(pvalues_0.01))
pfinal_0.01<-pchisq(test.stat_0.01, df=2*length(pvalues_0.01), lower.tail=F)
pvalues_0.1 <- c(pvalue_0.1_1,pvalue_0.1_2,pvalue_0.1_3,pvalue_0.1_4,pvalue_0.1_5,pvalue_0.1_6,pvalue_0.1_7, pvalue_0.1_8, pvalue_0.1_9, pvalue_0.1_13)
test.stat_0.1 <- -2*sum(log(pvalues_0.1))
pfinal_0.1<-pchisq(test.stat_0.1, df=2*length(pvalues_0.1), lower.tail=F)
pvalues_1.0 <- c(pvalue_1.0_1,pvalue_1.0_2,pvalue_1.0_3,pvalue_1.0_4,pvalue_1.0_5,pvalue_1.0_6,pvalue_1.0_7, pvalue_1.0_8, pvalue_1.0_9, pvalue_1.0_10, pvalue_1.0_11, pvalue_1.0_12, pvalue_1.0_13)
test.stat_1.0 <- -2*sum(log(pvalues_1.0))
pfinal_1.0<-pchisq(test.stat_1.0, df=2*length(pvalues_1.0), lower.tail=F)
pvalues_2.0 <- c(pvalue_2.0_3,pvalue_2.0_4,pvalue_2.0_5,pvalue_2.0_6,pvalue_2.0_7, pvalue_2.0_8, pvalue_2.0_9, pvalue_2.0_10, pvalue_2.0_11, pvalue_2.0_12, pvalue_2.0_13)
test.stat_2.0 <- -2*sum(log(pvalues_2.0))
pfinal_2.0<-pchisq(test.stat_2.0, df=2*length(pvalues_2.0), lower.tail=F)
pvalues_3.0 <- c(pvalue_3.0_10, pvalue_3.0_11, pvalue_3.0_12)
test.stat_3.0 <- -2*sum(log(pvalues_3.0))
pfinal_3.0<-pchisq(test.stat_3.0, df=2*length(pvalues_3.0), lower.tail=F)
pvalues_5.0 <- c(pvalue_5.0_8, pvalue_5.0_9, pvalue_5.0_10, pvalue_5.0_11, pvalue_5.0_12, pvalue_5.0_13)
test.stat_5.0 <- -2*sum(log(pvalues_5.0))
pfinal_5.0<-pchisq(test.stat_5.0, df=2*length(pvalues_5.0), lower.tail=F)
# Generate Table of final results
pvalues_0.01_table <- c(pvalue_0.01_1,pvalue_0.01_2,pvalue_0.01_3,pvalue_0.01_4,pvalue_0.01_5,pvalue_0.01_6,pvalue_0.01_7, NA, NA, NA, NA, NA, NA)
p0.01<-scales::pvalue(pvalues_0.01_table)
p0.01 <- p0.01 %>% replace_na("--")
pvalues_0.1_table <- c(pvalue_0.1_1,pvalue_0.1_2,pvalue_0.1_3,pvalue_0.1_4,pvalue_0.1_5,pvalue_0.1_6,pvalue_0.1_7, pvalue_0.1_8, pvalue_0.1_9,NA, NA, NA, pvalue_0.1_13)
p0.1<-scales::pvalue(pvalues_0.1_table)
p0.1<- p0.1 %>% replace_na("--")
pvalues_1.0_table <- c(pvalue_1.0_1,pvalue_1.0_2,pvalue_1.0_3,pvalue_1.0_4,pvalue_1.0_5,pvalue_1.0_6,pvalue_1.0_7, pvalue_1.0_8, pvalue_1.0_9, pvalue_1.0_10, pvalue_1.0_11, pvalue_1.0_12, pvalue_1.0_13)
p1.0<-scales::pvalue(pvalues_1.0_table)
pvalues_2.0_table <- c(NA, NA, pvalue_2.0_3, pvalue_2.0_4, pvalue_2.0_5, pvalue_2.0_6, pvalue_2.0_7, pvalue_2.0_8, pvalue_2.0_9, pvalue_2.0_10, pvalue_2.0_11, pvalue_2.0_12, pvalue_2.0_13)
p2.0<-scales::pvalue(pvalues_2.0_table)
p2.0<-p2.0 %>% replace_na("--")
pvalues_3.0_table <- c(NA, NA, NA, NA, NA, NA, NA, NA, NA, pvalue_3.0_10, pvalue_3.0_11, pvalue_3.0_12, NA)
p3.0<-scales::pvalue(pvalues_3.0_table)
p3.0<-p3.0 %>% replace_na("--")
pvalues_5.0_table <- c(NA, NA, NA, NA, NA, NA, NA, pvalue_5.0_8, pvalue_5.0_9, pvalue_5.0_10, pvalue_5.0_11, pvalue_5.0_12, pvalue_5.0_13)
p5.0<-scales::pvalue(pvalues_5.0_table)
p5.0<-p5.0 %>% replace_na("--")
pvalue_table <- data.frame("0.01" = p0.01, "0.1" = p0.1, "1.0" = p1.0, "2.0" = p2.0, "3.0" = p3.0, "5.0" = p5.0 )
colnames(pvalue_table) <- c(0.01, 0.1, 1.0, 2.0, 3.0, 5.0)
pfisher<- c(pfinal_0.01, pfinal_0.1, pfinal_1.0, pfinal_2.0, pfinal_3.0, pfinal_5.0)
pfisher<-scales::pvalue(pfisher)
pvalue_table[nrow(pvalue_table) + 1,]<-pfisher
row.names(pvalue_table) <- c("Block 1 - August 2nd, 2018", "Block 2 - August 6th, 2018", "Block 3 - June 14, 2019", "Block 4 - June 17, 2019", "Block 5 - June 19, 2019", "Block 6 - June 25, 2019", "Block 7 - June 27, 2019", "Block 8 - January 30, 2020", "Block 9 - February 13th, 2020", "Block 10 - September 27th, 2020", "Block 11 - October 6th, 2020", "Block 12 - November 1st, 2020", "Block 13 - June 16th, 2021", " Final P-value ")
pvalue_table %>%
kbl(caption = "<center><strong>P-values for log-rank tests on individual doses within each experiment<center><strong>") %>%
kable_classic(full_width = F, html_font = "Cambria", position = "center")%>%
add_header_above(c(" "=1, "Infectious Dose (Abs600nm)" = 6)) %>%
pack_rows("Individual Blocks", 1, 13) %>%
pack_rows("Fisher's Combined", 14, 14)
|
Infectious Dose (Abs600nm)
|
||||||
|---|---|---|---|---|---|---|
| 0.01 | 0.1 | 1 | 2 | 3 | 5 | |
| Individual Blocks | ||||||
| Block 1 - August 2nd, 2018 | 0.074 | 0.127 | 0.021 | – | – | – |
| Block 2 - August 6th, 2018 | 0.018 | 0.001 | 0.027 | – | – | – |
| Block 3 - June 14, 2019 | 0.065 | <0.001 | <0.001 | <0.001 | – | – |
| Block 4 - June 17, 2019 | 0.185 | 0.508 | <0.001 | <0.001 | – | – |
| Block 5 - June 19, 2019 | 0.205 | <0.001 | 0.003 | <0.001 | – | – |
| Block 6 - June 25, 2019 | 0.006 | <0.001 | <0.001 | <0.001 | – | – |
| Block 7 - June 27, 2019 | 0.650 | 0.008 | <0.001 | <0.001 | – | – |
| Block 8 - January 30, 2020 | – | <0.001 | <0.001 | <0.001 | – | <0.001 |
| Block 9 - February 13th, 2020 | – | <0.001 | <0.001 | <0.001 | – | <0.001 |
| Block 10 - September 27th, 2020 | – | – | <0.001 | <0.001 | <0.001 | <0.001 |
| Block 11 - October 6th, 2020 | – | – | <0.001 | <0.001 | <0.001 | <0.001 |
| Block 12 - November 1st, 2020 | – | – | <0.001 | <0.001 | <0.001 | <0.001 |
| Block 13 - June 16th, 2021 | – | <0.001 | <0.001 | <0.001 | – | <0.001 |
| Fisher’s Combined | ||||||
| Final P-value | <0.001 | <0.001 | <0.001 | <0.001 | <0.001 | <0.001 |
#KAPLAN MEIER
#table(Ef_Pr$primary, Ef_Pr$secondary_dose) #tells you the number of flies in each condition
#png(file="E.fae-P.rett.png", width=1200, heigh=900, res=130) #determines where the graph will be saved
par(mar=c(4,5,1,1)) #sets the margins around the graph so there is room for the labels
#The graph
plot(survfit(Surv(Ef_Pr$day, Ef_Pr$status) ~ Ef_Pr$primary + Ef_Pr$secondary_dose),
col=c("#B5E8FC", "#09B6FC", "#055df5", "#0600a8"), #line color, cycles through
lty=c(3, 3, 3, 3, 1, 1, 1, 1), #line type
lwd=5, #line width
yaxt='n', #hide the y axis so I can customize it
xaxt='n', #hide the x axis so I can customize itS
ylab="proportion alive", #label for the y-axis
xlab="days post-secondary infection", #label for the x-axis
cex.lab=2.5, #font size of axis labels
xlim=c(0,7), #determines range of the x-axis
xaxs='i') #makes the data flush with the axis
#Y axis
axis(2, #left side
las=2, #perpendicular to the axis
cex.axis=1.5, #font size
lwd=5) #line width
#X axis
axis(1, #bottom of the graph
cex.axis=1.5, #font size
lwd=5) #line width
box(lwd=5)
#The legend
legend("topright", #location
bty="n", #no box around the legend
c("PBS-P.rett 0.001", "PBS-P.rett 0.01", "PBS-P.rett 0.1", 'PBS-P.rett 1.0', "E.fae-P.rett 0.001", "E.fae-P.rett 0.01", "E.fae-P.rett 0.1", "E.fae-P.rett 1.0"), #labels in the legend
col=c("#B5E8FC", "#09B6FC", "#055df5", "#0600a8"), #colors for the legend lines, cycles through
lty=c(3, 3, 3, 3, 1, 1, 1, 1), #line type for legend lines
lwd=5, #line width for legend lines
cex=1) #font size
Ef_Pr_0.001<-subset(Ef_Pr,secondary_dose=="0.001")
#table(Ef_Pr_0.001$primary)
#THE GRAPH
#png(file="Figure2A.png",width=1200,height=900,res=130)
par(mar=c(6,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Ef_Pr_0.001), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=101)", "chronic & secondary infection (n=110)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.8)
#dev.off()
Ef_Pr_0.01<-subset(Ef_Pr,secondary_dose=="0.01")
#table(Ef_Pr_0.01$primary)
#THE GRAPH
#png(file="Figure2B.png",width=1200,height=900,res=130)
par(mar=c(6,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Ef_Pr_0.01), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=148)", "chronic & secondary infection (n=119)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.8)
#dev.off()
Ef_Pr_0.1<-subset(Ef_Pr,secondary_dose=="0.1")
#table(Ef_Pr_0.1$primary)
#THE GRAPH
#png(file="Figure2C.png",width=1200,height=900,res=130)
par(mar=c(6,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Ef_Pr_0.1), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=160)", "chronic & secondary infection (n=115)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.8)
#dev.off()
Ef_Pr_1.0<-subset(Ef_Pr,secondary_dose=="1")
#table(Ef_Pr_1.0$primary)
#THE GRAPH
#png(file="Figure2D.png",width=1200,height=900,res=130)
par(mar=c(6,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Ef_Pr_1.0), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=157)", "chronic & secondary infection (n=125)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.8)
#dev.off()
g <- ggdraw() +
draw_image("Figure2A.png")
h <- ggdraw() +
draw_image("Figure2B.png")
i <- ggdraw() +
draw_image("Figure2C.png")
j <- ggdraw() +
draw_image("Figure2D.png")
multiplot <- g + h + i + j
y <- multiplot +
plot_annotation(tag_levels = 'A')
ggsave("Fig.2.ChronicEf_Pr.pdf", y, device = "pdf")
y
Ef_Pr_block1 <- subset(Ef_Pr, Date == "6/11/2019")
table(Ef_Pr_block1 $primary, Ef_Pr_block1 $secondary_dose)
##
## 0.01 0.1 1
## 01_PBS 20 19 20
## 02_Ef 16 16 17
#Ef_Pr_0.001_1<-subset(Ef_Pr_block1, secondary_dose=="2")
#table(Ef_Pr_0.001_1$primary)
Ef_Pr_0.01_1<-subset(Ef_Pr_block1, secondary_dose=="0.01")
#table(Ef_Pr_0.01_1$primary)
Ef_Pr_0.1_1<-subset(Ef_Pr_block1, secondary_dose=="0.1")
#table(Ef_Pr_0.1_1$primary)
Ef_Pr_1.0_1<-subset(Ef_Pr_block1, secondary_dose=="1")
#table(Ef_Pr_1.0_1$primary)
surv_Ef_Pr_0.01_1<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.01_1)
surv_Ef_Pr_0.01_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.01_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 20 19 12.8 3.03 8.21
## primary=02_Ef 16 16 22.2 1.74 8.21
##
## Chisq= 8.2 on 1 degrees of freedom, p= 0.004
surv_Ef_Pr_0.1_1<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.1_1)
surv_Ef_Pr_0.1_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.1_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 19 19 14.2 1.66 5.9
## primary=02_Ef 16 15 19.8 1.18 5.9
##
## Chisq= 5.9 on 1 degrees of freedom, p= 0.02
surv_Ef_Pr_1.0_1<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_1.0_1)
surv_Ef_Pr_1.0_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_1.0_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 20 20 17.8 0.262 5.13
## primary=02_Ef 17 17 19.2 0.244 5.13
##
## Chisq= 5.1 on 1 degrees of freedom, p= 0.02
par(mar=c(4,5,1,1))
plot(survfit(Surv(Ef_Pr_block1$day, Ef_Pr_block1$status) ~ Ef_Pr_block1$primary + Ef_Pr_block1$secondary_dose), col=c("#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.01", "PBS-P.rett 0.1", 'PBS-P.rett 1.0', "E.fae-P.rett 0.01", "E.fae-P.rett 0.1", "E.fae-P.rett 1.0"),col=c("#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 1, 1, 1),lwd=5, cex=1)
Ef_Pr_block2 <- subset(Ef_Pr, Date == "6/22/2018")
table(Ef_Pr_block2 $primary, Ef_Pr_block2 $secondary_dose)
##
## 0.001 0.01 0.1 1
## 01_PBS 28 18 29 30
## 02_Ef 33 18 19 29
Ef_Pr_0.001_2<-subset(Ef_Pr_block2, secondary_dose=="0.001")
table(Ef_Pr_0.001_2$primary)
##
## 01_PBS 02_Ef
## 28 33
Ef_Pr_0.01_2<-subset(Ef_Pr_block2, secondary_dose=="0.01")
#table(Ef_Pr_0.01_2$primary)
Ef_Pr_0.1_2<-subset(Ef_Pr_block2, secondary_dose=="0.1")
#table(Ef_Pr_0.1_2$primary)
Ef_Pr_1.0_2<-subset(Ef_Pr_block2, secondary_dose=="1")
#table(Ef_Pr_1.0_2$primary)
surv_Ef_Pr_0.001_2<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.001_2)
surv_Ef_Pr_0.001_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.001_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 28 13 10.9 0.406 1.01
## primary=02_Ef 33 11 13.1 0.337 1.01
##
## Chisq= 1 on 1 degrees of freedom, p= 0.3
surv_Ef_Pr_0.01_2<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.01_2)
surv_Ef_Pr_0.01_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.01_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 18 13 8.44 2.46 6.03
## primary=02_Ef 18 6 10.56 1.97 6.03
##
## Chisq= 6 on 1 degrees of freedom, p= 0.01
surv_Ef_Pr_0.1_2<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.1_2)
surv_Ef_Pr_0.1_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.1_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 29 23 17.1 2.01 7.88
## primary=02_Ef 19 8 13.9 2.48 7.88
##
## Chisq= 7.9 on 1 degrees of freedom, p= 0.005
surv_Ef_Pr_1.0_2<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_1.0_2)
surv_Ef_Pr_1.0_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_1.0_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 30 30 27.5 0.235 5.56
## primary=02_Ef 29 26 28.5 0.226 5.56
##
## Chisq= 5.6 on 1 degrees of freedom, p= 0.02
par(mar=c(4,5,1,1))
plot(survfit(Surv(Ef_Pr_block2$day, Ef_Pr_block2$status) ~ Ef_Pr_block2$primary + Ef_Pr_block2$secondary_dose), col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 3, 1, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.001", "PBS-P.rett 0.01", "PBS-P.rett 0.1", 'PBS-P.rett 1.0', "E.fae-P.rett 0.001", "E.fae-P.rett 0.01", "E.fae-P.rett 0.1", "E.fae-P.rett 1.0"),col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 3, 1, 1, 1, 1),lwd=5, cex=1)
Ef_Pr_block3 <- subset(Ef_Pr, Date == "6/25/2018")
table(Ef_Pr_block3 $primary, Ef_Pr_block3 $secondary_dose)
##
## 0.001 0.01 0.1 1
## 01_PBS 25 27 29 37
## 02_Ef 23 15 25 19
Ef_Pr_0.001_3<-subset(Ef_Pr_block3, secondary_dose=="0.001")
table(Ef_Pr_0.001_3$primary)
##
## 01_PBS 02_Ef
## 25 23
Ef_Pr_0.01_3<-subset(Ef_Pr_block3, secondary_dose=="0.01")
#table(Ef_Pr_0.01_3$primary)
Ef_Pr_0.1_3<-subset(Ef_Pr_block3, secondary_dose=="0.1")
#table(Ef_Pr_0.1_3$primary)
Ef_Pr_1.0_3<-subset(Ef_Pr_block3, secondary_dose=="1")
#table(Ef_Pr_1.0_3$primary)
surv_Ef_Pr_0.001_3<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.001_3)
surv_Ef_Pr_0.001_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.001_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 25 18 14.4 0.912 3.13
## primary=02_Ef 23 11 14.6 0.897 3.13
##
## Chisq= 3.1 on 1 degrees of freedom, p= 0.08
surv_Ef_Pr_0.01_3<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.01_3)
surv_Ef_Pr_0.01_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.01_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 27 21 15.8 1.68 7.31
## primary=02_Ef 15 5 10.2 2.62 7.31
##
## Chisq= 7.3 on 1 degrees of freedom, p= 0.007
surv_Ef_Pr_0.1_3<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.1_3)
surv_Ef_Pr_0.1_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.1_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 29 21 20.6 0.00806 0.0384
## primary=02_Ef 25 18 18.4 0.00902 0.0384
##
## Chisq= 0 on 1 degrees of freedom, p= 0.8
surv_Ef_Pr_1.0_3<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_1.0_3)
surv_Ef_Pr_1.0_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_1.0_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 37 36 35.7 0.00290 0.235
## primary=02_Ef 19 18 18.3 0.00564 0.235
##
## Chisq= 0.2 on 1 degrees of freedom, p= 0.6
par(mar=c(4,5,1,1))
plot(survfit(Surv(Ef_Pr_block3$day, Ef_Pr_block3$status) ~ Ef_Pr_block3$primary + Ef_Pr_block3$secondary_dose), col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 3, 1, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.001", "PBS-P.rett 0.01", "PBS-P.rett 0.1", 'PBS-P.rett 1.0', "E.fae-P.rett 0.001", "E.fae-P.rett 0.01", "E.fae-P.rett 0.1", "E.fae-P.rett 1.0"),col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 3, 1, 1, 1, 1),lwd=5, cex=1)
Not sure why the data stops after 3 days, Frank did this one
Ef_Pr_block4 <- subset(Ef_Pr, Date == "7/19/2018")
table(Ef_Pr_block4 $primary, Ef_Pr_block4 $secondary_dose)
##
## 0.001 0.01 0.1 1
## 01_PBS 9 18 18 10
## 02_Ef 23 27 16 17
Ef_Pr_0.001_4<-subset(Ef_Pr_block4, secondary_dose=="0.001")
#table(Ef_Pr_0.001_4$primary)
Ef_Pr_0.01_4<-subset(Ef_Pr_block4, secondary_dose=="0.01")
#table(Ef_Pr_0.01_4$primary)
Ef_Pr_0.1_4<-subset(Ef_Pr_block4, secondary_dose=="0.1")
#table(Ef_Pr_0.1_4$primary)
Ef_Pr_1.0_4<-subset(Ef_Pr_block4, secondary_dose=="1")
#table(Ef_Pr_1.0_4$primary)
surv_Ef_Pr_0.001_4<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.001_4)
surv_Ef_Pr_0.001_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.001_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 9 3 2.56 0.0743 0.123
## primary=02_Ef 23 6 6.44 0.0296 0.123
##
## Chisq= 0.1 on 1 degrees of freedom, p= 0.7
surv_Ef_Pr_0.01_4<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.01_4)
surv_Ef_Pr_0.01_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.01_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 18 12 12 7.93e-05 0.00026
## primary=02_Ef 27 16 16 5.98e-05 0.00026
##
## Chisq= 0 on 1 degrees of freedom, p= 1
surv_Ef_Pr_0.1_4<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.1_4)
surv_Ef_Pr_0.1_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.1_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 18 9 8.31 0.0574 0.147
## primary=02_Ef 16 9 9.69 0.0492 0.147
##
## Chisq= 0.1 on 1 degrees of freedom, p= 0.7
surv_Ef_Pr_1.0_4<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_1.0_4)
surv_Ef_Pr_1.0_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_1.0_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 10 9 8.79 0.00516 0.0292
## primary=02_Ef 17 17 17.21 0.00263 0.0292
##
## Chisq= 0 on 1 degrees of freedom, p= 0.9
par(mar=c(4,5,1,1))
plot(survfit(Surv(Ef_Pr_block4$day, Ef_Pr_block4$status) ~ Ef_Pr_block4$primary + Ef_Pr_block4$secondary_dose), col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 3, 1, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.001", "PBS-P.rett 0.01", "PBS-P.rett 0.1", 'PBS-P.rett 1.0', "E.fae-P.rett 0.001", "E.fae-P.rett 0.01", "E.fae-P.rett 0.1", "E.fae-P.rett 1.0"),col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 3, 1, 1, 1, 1),lwd=5, cex=1)
Ef_Pr_block5 <- subset(Ef_Pr, Date == "6/3/2019")
table(Ef_Pr_block5 $primary, Ef_Pr_block5 $secondary_dose)
##
## 0.001 0.01 0.1 1
## 01_PBS 18 21 19 17
## 02_Ef 12 13 6 12
Ef_Pr_0.001_5<-subset(Ef_Pr_block5, secondary_dose=="0.001")
table(Ef_Pr_0.001_5$primary)
##
## 01_PBS 02_Ef
## 18 12
Ef_Pr_0.01_5<-subset(Ef_Pr_block5, secondary_dose=="0.01")
#table(Ef_Pr_0.01_5$primary)
Ef_Pr_0.1_5<-subset(Ef_Pr_block5, secondary_dose=="0.1")
#table(Ef_Pr_0.1_5$primary)
Ef_Pr_1.0_5<-subset(Ef_Pr_block5, secondary_dose=="1")
#table(Ef_Pr_1.0_5$primary)
surv_Ef_Pr_0.001_5<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.001_5)
surv_Ef_Pr_0.001_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.001_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 18 13 10.65 0.518 1.39
## primary=02_Ef 12 6 8.35 0.661 1.39
##
## Chisq= 1.4 on 1 degrees of freedom, p= 0.2
surv_Ef_Pr_0.01_5<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.01_5)
surv_Ef_Pr_0.01_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.01_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 21 16 12.9 0.759 2.31
## primary=02_Ef 13 7 10.1 0.965 2.31
##
## Chisq= 2.3 on 1 degrees of freedom, p= 0.1
surv_Ef_Pr_0.1_5<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.1_5)
surv_Ef_Pr_0.1_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.1_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 19 15 14.05 0.0636 0.502
## primary=02_Ef 6 4 4.95 0.1806 0.502
##
## Chisq= 0.5 on 1 degrees of freedom, p= 0.5
surv_Ef_Pr_1.0_5<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_1.0_5)
surv_Ef_Pr_1.0_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_1.0_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 17 17 14.7 0.375 6.35
## primary=02_Ef 12 12 14.3 0.383 6.35
##
## Chisq= 6.3 on 1 degrees of freedom, p= 0.01
par(mar=c(4,5,1,1))
plot(survfit(Surv(Ef_Pr_block5$day, Ef_Pr_block5$status) ~ Ef_Pr_block5$primary + Ef_Pr_block5$secondary_dose), col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 3, 1, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.001", "PBS-P.rett 0.01", "PBS-P.rett 0.1", 'PBS-P.rett 1.0', "E.fae-P.rett 0.001", "E.fae-P.rett 0.01", "E.fae-P.rett 0.1", "E.fae-P.rett 1.0"),col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 3, 1, 1, 1, 1),lwd=5, cex=1)
Ef_Pr_block6 <- subset(Ef_Pr, Date == "6/5/2019")
table(Ef_Pr_block6 $primary, Ef_Pr_block6 $secondary_dose)
##
## 0.001 0.01 0.1 1
## 01_PBS 21 21 22 20
## 02_Ef 19 16 16 17
Ef_Pr_0.001_6<-subset(Ef_Pr_block6, secondary_dose=="0.001")
table(Ef_Pr_0.001_6$primary)
##
## 01_PBS 02_Ef
## 21 19
Ef_Pr_0.01_6<-subset(Ef_Pr_block6, secondary_dose=="0.01")
#table(Ef_Pr_0.01_6$primary)
Ef_Pr_0.1_6<-subset(Ef_Pr_block6, secondary_dose=="0.1")
#table(Ef_Pr_0.1_6$primary)
Ef_Pr_1.0_6<-subset(Ef_Pr_block6, secondary_dose=="1")
#table(Ef_Pr_1.0_6$primary)
surv_Ef_Pr_0.001_6<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.001_6)
surv_Ef_Pr_0.001_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.001_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 21 13 8.18 2.84 6.64
## primary=02_Ef 19 4 8.82 2.63 6.64
##
## Chisq= 6.6 on 1 degrees of freedom, p= 0.01
surv_Ef_Pr_0.01_6<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.01_6)
surv_Ef_Pr_0.01_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.01_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 21 10 7.63 0.738 1.76
## primary=02_Ef 16 4 6.37 0.883 1.76
##
## Chisq= 1.8 on 1 degrees of freedom, p= 0.2
surv_Ef_Pr_0.1_6<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.1_6)
surv_Ef_Pr_0.1_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.1_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 22 10 10.96 0.0834 0.246
## primary=02_Ef 16 9 8.04 0.1137 0.246
##
## Chisq= 0.2 on 1 degrees of freedom, p= 0.6
surv_Ef_Pr_1.0_6<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_1.0_6)
surv_Ef_Pr_1.0_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_1.0_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 20 17 16.7 0.00591 0.026
## primary=02_Ef 17 16 16.3 0.00604 0.026
##
## Chisq= 0 on 1 degrees of freedom, p= 0.9
par(mar=c(4,5,1,1))
plot(survfit(Surv(Ef_Pr_block6$day, Ef_Pr_block6$status) ~ Ef_Pr_block6$primary + Ef_Pr_block6$secondary_dose), col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 3, 1, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.001", "PBS-P.rett 0.01", "PBS-P.rett 0.1", 'PBS-P.rett 1.0', "E.fae-P.rett 0.001", "E.fae-P.rett 0.01", "E.fae-P.rett 0.1", "E.fae-P.rett 1.0"),col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 3, 1, 1, 1, 1),lwd=5, cex=1)
Ef_Pr_block7 <- subset(Ef_Pr, Date == "7/1/2019")
table(Ef_Pr_block7 $primary, Ef_Pr_block7 $secondary_dose)
##
## 0.01 0.1 1
## 01_PBS 23 24 23
## 02_Ef 14 17 14
#Ef_Pr_0.001_7<-subset(Ef_Pr_block3, secondary_dose=="0.001")
#table(Ef_Pr_0.001_7$primary)
Ef_Pr_0.01_7<-subset(Ef_Pr_block3, secondary_dose=="0.01")
#table(Ef_Pr_0.01_7$primary)
Ef_Pr_0.1_7<-subset(Ef_Pr_block3, secondary_dose=="0.1")
#table(Ef_Pr_0.1_7$primary)
Ef_Pr_1.0_7<-subset(Ef_Pr_block3, secondary_dose=="1")
#table(Ef_Pr_1.0_7$primary)
#surv_Ef_Pr_0.001_7<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.001_7)
#surv_Ef_Pr_0.001_7
surv_Ef_Pr_0.01_7<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.01_7)
surv_Ef_Pr_0.01_7
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.01_7)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 27 21 15.8 1.68 7.31
## primary=02_Ef 15 5 10.2 2.62 7.31
##
## Chisq= 7.3 on 1 degrees of freedom, p= 0.007
surv_Ef_Pr_0.1_7<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.1_7)
surv_Ef_Pr_0.1_7
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.1_7)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 29 21 20.6 0.00806 0.0384
## primary=02_Ef 25 18 18.4 0.00902 0.0384
##
## Chisq= 0 on 1 degrees of freedom, p= 0.8
surv_Ef_Pr_1.0_7<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_1.0_7)
surv_Ef_Pr_1.0_7
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_1.0_7)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 37 36 35.7 0.00290 0.235
## primary=02_Ef 19 18 18.3 0.00564 0.235
##
## Chisq= 0.2 on 1 degrees of freedom, p= 0.6
par(mar=c(4,5,1,1))
plot(survfit(Surv(Ef_Pr_block7$day, Ef_Pr_block7$status) ~ Ef_Pr_block7$primary + Ef_Pr_block7$secondary_dose), col=c("#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.01", "PBS-P.rett 0.1", 'PBS-P.rett 1.0', "E.fae-P.rett 0.01", "E.fae-P.rett 0.1", "E.fae-P.rett 1.0"),col=c("#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 1, 1, 1),lwd=5, cex=1)
# Grab p-values from each individual logrank test (results shown above in each block)
pvalue_Ef_Pr_0.01_1<-pchisq(surv_Ef_Pr_0.01_1$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.1_1<-pchisq(surv_Ef_Pr_0.1_1$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_1.0_1<-pchisq(surv_Ef_Pr_1.0_1$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.001_2<-pchisq(surv_Ef_Pr_0.001_2$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.01_2<-pchisq(surv_Ef_Pr_0.01_2$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.1_2<-pchisq(surv_Ef_Pr_0.1_2$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_1.0_2<-pchisq(surv_Ef_Pr_1.0_2$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.001_3<-pchisq(surv_Ef_Pr_0.001_3$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.01_3<-pchisq(surv_Ef_Pr_0.01_3$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.1_3<-pchisq(surv_Ef_Pr_0.1_3$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_1.0_3<-pchisq(surv_Ef_Pr_1.0_3$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.001_4<-pchisq(surv_Ef_Pr_0.001_4$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.01_4<-pchisq(surv_Ef_Pr_0.01_4$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.1_4<-pchisq(surv_Ef_Pr_0.1_4$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_1.0_4<-pchisq(surv_Ef_Pr_1.0_4$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.001_5<-pchisq(surv_Ef_Pr_0.001_5$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.01_5<-pchisq(surv_Ef_Pr_0.01_5$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.1_5<-pchisq(surv_Ef_Pr_0.1_5$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_1.0_5<-pchisq(surv_Ef_Pr_1.0_5$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.001_6<-pchisq(surv_Ef_Pr_0.001_6$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.01_6<-pchisq(surv_Ef_Pr_0.01_6$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.1_6<-pchisq(surv_Ef_Pr_0.1_6$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_1.0_6<-pchisq(surv_Ef_Pr_1.0_6$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.01_7<-pchisq(surv_Ef_Pr_0.01_7$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.1_7<-pchisq(surv_Ef_Pr_0.1_7$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_1.0_7<-pchisq(surv_Ef_Pr_1.0_7$chisq, df=1, lower.tail=F)
# Use all p-values from a given condition to get Fisher's combined test p-value# Use all p-values from a given condition to get Fisher's combined test p-value
pvalues_Ef_Pr_0.001 <- c(pvalue_Ef_Pr_0.001_2,pvalue_Ef_Pr_0.001_3,pvalue_Ef_Pr_0.001_4,pvalue_Ef_Pr_0.001_5,pvalue_Ef_Pr_0.001_6)
test.stat_Ef_Pr_0.001 <- -2*sum(log(pvalues_Ef_Pr_0.001))
pfinal_Ef_Pr_0.001<-pchisq(test.stat_Ef_Pr_0.001, df=2*length(pvalues_Ef_Pr_0.001), lower.tail=F)
pvalues_Ef_Pr_0.01 <- c(pvalue_Ef_Pr_0.01_1,pvalue_Ef_Pr_0.01_2,pvalue_Ef_Pr_0.01_3,pvalue_Ef_Pr_0.01_4,pvalue_Ef_Pr_0.01_5,pvalue_Ef_Pr_0.01_6, pvalue_Ef_Pr_0.01_7)
test.stat_Ef_Pr_0.01 <- -2*sum(log(pvalues_Ef_Pr_0.01))
pfinal_Ef_Pr_0.01<-pchisq(test.stat_Ef_Pr_0.01, df=2*length(pvalues_Ef_Pr_0.01), lower.tail=F)
pvalues_Ef_Pr_0.1 <- c(pvalue_Ef_Pr_0.1_1,pvalue_Ef_Pr_0.1_2,pvalue_Ef_Pr_0.1_3,pvalue_Ef_Pr_0.1_4,pvalue_Ef_Pr_0.1_5,pvalue_Ef_Pr_0.1_6,pvalue_Ef_Pr_0.1_7)
test.stat_Ef_Pr_0.1 <- -2*sum(log(pvalues_Ef_Pr_0.1))
pfinal_Ef_Pr_0.1<-pchisq(test.stat_Ef_Pr_0.1, df=2*length(pvalues_Ef_Pr_0.1), lower.tail=F)
pvalues_Ef_Pr_1.0 <- c(pvalue_Ef_Pr_1.0_1,pvalue_Ef_Pr_1.0_2,pvalue_Ef_Pr_1.0_3,pvalue_Ef_Pr_1.0_4,pvalue_Ef_Pr_1.0_5,pvalue_Ef_Pr_1.0_6,pvalue_Ef_Pr_1.0_7)
test.stat_Ef_Pr_1.0 <- -2*sum(log(pvalues_Ef_Pr_1.0))
pfinal_Ef_Pr_1.0<-pchisq(test.stat_Ef_Pr_1.0, df=2*length(pvalues_Ef_Pr_1.0), lower.tail=F)
#### Generate Table of final results
pvalues_Ef_Pr_0.001_table <- c(NA,pvalue_Ef_Pr_0.001_2,pvalue_Ef_Pr_0.001_3,pvalue_Ef_Pr_0.001_4,pvalue_Ef_Pr_0.001_5,pvalue_Ef_Pr_0.001_6,NA)
p_Ef_Pr_0.001<-scales::pvalue(pvalues_Ef_Pr_0.001_table)
p_Ef_Pr_0.001 <- p_Ef_Pr_0.001 %>% replace_na("--")
pvalues_Ef_Pr_0.01_table <- c(pvalue_Ef_Pr_0.01_1,pvalue_Ef_Pr_0.01_2,pvalue_Ef_Pr_0.01_3,pvalue_Ef_Pr_0.01_4,pvalue_Ef_Pr_0.01_5,pvalue_Ef_Pr_0.01_6,pvalue_Ef_Pr_0.01_7)
p_Ef_Pr_0.01<-scales::pvalue(pvalues_Ef_Pr_0.01_table)
pvalues_Ef_Pr_0.1_table <- c(pvalue_Ef_Pr_0.1_1,pvalue_Ef_Pr_0.1_2,pvalue_Ef_Pr_0.1_3,pvalue_Ef_Pr_0.1_4,pvalue_Ef_Pr_0.1_5,pvalue_Ef_Pr_0.1_6,pvalue_Ef_Pr_0.1_7)
p_Ef_Pr_0.1<-scales::pvalue(pvalues_Ef_Pr_0.1_table)
pvalues_Ef_Pr_1.0_table <- c(pvalue_Ef_Pr_1.0_1,pvalue_Ef_Pr_1.0_2,pvalue_Ef_Pr_1.0_3,pvalue_Ef_Pr_1.0_4,pvalue_Ef_Pr_1.0_5,pvalue_Ef_Pr_1.0_6,pvalue_Ef_Pr_1.0_7)
p_Ef_Pr_1.0<-scales::pvalue(pvalues_Ef_Pr_1.0_table)
pvalue_table_Ef_Pr <- data.frame("0.001" = p_Ef_Pr_0.001, "0.01" = p_Ef_Pr_0.01, "0.1" = p_Ef_Pr_0.1, "1.0" = p_Ef_Pr_1.0)
colnames(pvalue_table_Ef_Pr) <- c(0.001,0.01, 0.1, 1.0)
pfisher_Ef_Pr<- c(pfinal_Ef_Pr_0.001, pfinal_Ef_Pr_0.01, pfinal_Ef_Pr_0.1, pfinal_Ef_Pr_1.0)
pfisher_Ef_Pr<-scales::pvalue(pfisher_Ef_Pr)
pvalue_table_Ef_Pr[nrow(pvalue_table_Ef_Pr) + 1,]<-pfisher_Ef_Pr
row.names(pvalue_table_Ef_Pr) <- c("Block 1 - June 11th, 2018", "Block 2 - June 22nd, 2018", "Block 3 - June 25th, 2018", "Block 4 - July 19th, 2018", "Block 5 - June 3, 2019", "Block 6 - June 5, 2019", "Block 7 - July 1st, 2019", " Final P-value ")
pvalue_table_Ef_Pr %>%
kbl(caption = "<center><strong>P-values for log-rank tests on individual doses within each experiment<center><strong>") %>%
kable_classic(full_width = F, html_font = "Cambria", position = "center")%>%
add_header_above(c(" "=1, "Infectious Dose (Abs600nm)" = 4)) %>%
pack_rows("Individual Blocks", 1, 7) %>%
pack_rows("Fisher's Combined", 8, 8)
|
Infectious Dose (Abs600nm)
|
||||
|---|---|---|---|---|
| 0.001 | 0.01 | 0.1 | 1 | |
| Individual Blocks | ||||
| Block 1 - June 11th, 2018 | – | 0.004 | 0.015 | 0.023 |
| Block 2 - June 22nd, 2018 | 0.315 | 0.014 | 0.005 | 0.018 |
| Block 3 - June 25th, 2018 | 0.077 | 0.007 | 0.845 | 0.628 |
| Block 4 - July 19th, 2018 | 0.725 | 0.987 | 0.702 | 0.864 |
| Block 5 - June 3, 2019 | 0.239 | 0.129 | 0.479 | 0.012 |
| Block 6 - June 5, 2019 | 0.010 | 0.184 | 0.620 | 0.872 |
| Block 7 - July 1st, 2019 | – | 0.007 | 0.845 | 0.628 |
| Fisher’s Combined | ||||
| Final P-value | 0.028 | <0.001 | 0.064 | 0.020 |
Study design:
#png(file="Figure3A.png", width=3200, height=1800, res=300)
ggplot(Sm_bacterial_load, aes(x = as.character(secondary_dose), y = log(load_10hr), fill = primary)) +
scale_fill_manual(values=c("white", "#3399FF"),
name=" ",
breaks=c("01_PBS","02_Sm"),
labels=c("control", "chronically infected")) +
geom_boxplot() +
#ggtitle("S. marcescens chronic infection") +
theme_classic(base_size = 24) +
xlab("secondary dose (Abs600nm)") +
ylab("log(calculated bacterial load)") +
guides(fill=guide_legend(title=NULL)) +
geom_point(color = "black", size = 1.3, position=position_jitterdodge(jitter.width=0.2))
#dev.off()
S.m_resistance <- lm(log(load_10hr) ~ rep + secondary_dose + primary + secondary_dose:primary, data = Sm_bacterial_load)
S.m_resistance_anova <- anova(S.m_resistance)
S.m_resistance_anova
## Analysis of Variance Table
##
## Response: log(load_10hr)
## Df Sum Sq Mean Sq F value Pr(>F)
## rep 8 56.729 7.091 2.8793 0.01059 *
## secondary_dose 1 145.561 145.561 59.1044 6.543e-10 ***
## primary 1 94.077 94.077 38.1995 1.328e-07 ***
## secondary_dose:primary 1 2.651 2.651 1.0765 0.30467
## Residuals 48 118.213 2.463
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
shapiro.test(resid(S.m_resistance))
##
## Shapiro-Wilk normality test
##
## data: resid(S.m_resistance)
## W = 0.98595, p-value = 0.7192
#BACTERIAL LOAD
Ef_bacterial_load_10hr <- subset(Ef_bacterial_load, time == 10)
#png(file="Figure3B.png",width=3200, height=1800, res=300)
ggplot(Ef_bacterial_load_10hr, aes(x = as.character(secondary_dose), y = log(load), fill = primary)) +
scale_fill_manual(values=c("white", "#3399FF"),
name=" ",
breaks=c("01_PBS","02_Efae"),
labels=c("control", "chronically infected")) +
geom_boxplot() +
#ggtitle("E. faecalis chronic infection") +
theme_classic(base_size = 24) +
xlab("secondary dose (Abs600nm)") +
ylab("log(CFU/fly)") +
guides(fill=guide_legend(title=NULL)) +
geom_point(color = "black", size = 1.3, position=position_jitterdodge(jitter.width=0.2))
#dev.off()
E.fae_resistance <- lm(log(load) ~ replicate + secondary_dose + primary + secondary_dose:primary, data = Ef_bacterial_load_10hr)
E.fae_resistance_anova <- anova(E.fae_resistance)
E.fae_resistance_anova
## Analysis of Variance Table
##
## Response: log(load)
## Df Sum Sq Mean Sq F value Pr(>F)
## replicate 5 231.74 46.35 15.403 3.504e-13 ***
## secondary_dose 1 799.17 799.17 265.589 < 2.2e-16 ***
## primary 1 35.09 35.09 11.660 0.0007463 ***
## secondary_dose:primary 1 0.23 0.23 0.075 0.7843830
## Residuals 247 743.23 3.01
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
shapiro.test(resid(E.fae_resistance))
##
## Shapiro-Wilk normality test
##
## data: resid(E.fae_resistance)
## W = 0.99623, p-value = 0.7984
k <- ggdraw() +
ggtitle('P. rettgeri') +
draw_image("Figure3A.png")
l <- ggdraw() +
draw_image("Figure3B.png")
m <- ggdraw() +
draw_image("Figure3C.png")
n <- ggdraw() +
draw_image("Figure3D.png")
multiplot3 <- k + m + l + n
x<-multiplot3
#plot_annotation(tag_levels = 'A')
ggsave("Fig.3.ResistanceTolerance.pdf", x, device = "pdf")
x
Study design:
#table(Sm_Ps$primary, Sm_Ps$secondary_dose)
#png(file="Sm_Ps.png", width = 1700, height = 1500, res = 300)
#par(mar=c(4,4,0.1,0.1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(Sm_Ps$hour, Sm_Ps$status) ~ Sm_Ps$primary + Sm_Ps$secondary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1), lwd=2.6, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Hours Post-Secondary Infection", cex.lab=1.3, xlim=c(10,48), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-PBS", "PBS-0.001", "PBS-0.01", "PBS-0.1", "PBS-1.0", "S.m-PBS", "S.m-0.001", "S.m-0.01", "S.m-0.1", "S.m-1.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1),lwd=4, cex=1)
Sm_Ps_PBS<-subset(Sm_Ps,secondary_dose=="0")
#table(Sm_Ps_PBS$primary)
#THE GRAPH
#png(file="AFigure1B.png",width=1200,height=900,res=130)
par(mar=c(4,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Sm_Ps_PBS), col=c("gray"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("bottomleft",inset=0, bty="n" , c("double injection control (n=59)", "chronic infection only (n=59)"),col=c("gray"),lty=c(3, 1),lwd=5, cex=1)
#dev.off()
Sm_Ps_0.001<-subset(Sm_Ps,secondary_dose=="0.001")
#table(Sm_Ps_0.001$primary)
#THE GRAPH
#png(file="Figure4A.png",width=1200,height=900,res=130)
par(mar=c(6,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(hour, status) ~ primary + secondary_dose, data=Sm_Ps_0.001), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="hours post-secondary infection", cex.lab=2, xlim=c(0,72), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("bottomleft",inset=0, bty="n" , c("secondary infection only (n=82)", "chronic & secondary infection (n=66)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.6)
#dev.off()
Sm_Ps_0.01<-subset(Sm_Ps,secondary_dose=="0.01")
#table(Sm_Ps_0.01$primary)
#THE GRAPH
#png(file="Figure4B.png",width=1200,height=900,res=130)
par(mar=c(6,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(hour, status) ~ primary + secondary_dose, data=Sm_Ps_0.01), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="hours post-secondary infection", cex.lab=2, xlim=c(0,72), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("bottomleft",inset=0, bty="n" , c("secondary infection only (n=87)", "chronic & secondary infection (n=80)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.6)
#dev.off()
Sm_Ps_0.1<-subset(Sm_Ps,secondary_dose=="0.1")
#table(Sm_Ps_0.1$primary)
#THE GRAPH
#png(file="Figure4C.png",width=1200,height=900,res=130)
par(mar=c(6,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(hour, status) ~ primary + secondary_dose, data=Sm_Ps_0.1), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="hours post-secondary infection", cex.lab=2, xlim=c(0,72), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=80)", "chronic & secondary infection (n=81)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.6)
#dev.off()
Sm_Ps_1.0<-subset(Sm_Ps,secondary_dose=="1")
#table(Sm_Ps_1.0$primary)
#png(file="Figure4D.png",width=1200,height=900,res=130)
par(mar=c(6,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(hour, status) ~ primary + secondary_dose, data=Sm_Ps_1.0), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="hours post-secondary infection", cex.lab=2, xlim=c(0,72), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=75)", "chronic & secondary infection (n=81)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.6)
#dev.off()
m <- ggdraw() +
draw_image("Figure4A.png")
n <- ggdraw() +
draw_image("Figure4B.png")
o <- ggdraw() +
ggtitle('P. rettgeri') +
draw_image("Figure4C.png")
p <- ggdraw() +
draw_image("Figure4D.png")
multiplot3 <- m + n + o + p
w<-multiplot3 +
plot_annotation(tag_levels = 'A')
ggsave("Fig.4.Psneebia.pdf", w, device = "pdf")
w
Sm_Ps_block1 <- subset(Sm_Ps, date == "10/1/2020")
table(Sm_Ps_block1$primary, Sm_Ps_block1$secondary_dose)
##
## 0 0.001 0.01 0.1 1
## PBS 19 30 29 20 30
## S.m 19 16 24 29 25
#subset by secondary dose
Sm_Ps_0.001_1<-subset(Sm_Ps_block1, secondary_dose=="0.001")
Sm_Ps_0.01_1<-subset(Sm_Ps_block1, secondary_dose=="0.01")
Sm_Ps_0.1_1<-subset(Sm_Ps_block1, secondary_dose=="0.1")
Sm_Ps_1.0_1<-subset(Sm_Ps_block1, secondary_dose=="1")
surv_Sm_Ps_0.001_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_0.001_1)
surv_Sm_Ps_0.001_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_0.001_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30 21 18.3 0.388 1.13
## primary=S.m 16 10 12.7 0.562 1.13
##
## Chisq= 1.1 on 1 degrees of freedom, p= 0.3
surv_Sm_Ps_0.01_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_0.01_1)
surv_Sm_Ps_0.01_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_0.01_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 29 29 14.7 13.87 29.7
## primary=S.m 24 23 37.3 5.47 29.7
##
## Chisq= 29.7 on 1 degrees of freedom, p= 5e-08
surv_Sm_Ps_0.1_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_0.1_1)
surv_Sm_Ps_0.1_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_0.1_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 20 20 6.8 25.65 45.2
## primary=S.m 29 29 42.2 4.13 45.2
##
## Chisq= 45.2 on 1 degrees of freedom, p= 2e-11
surv_Sm_Ps_1.0_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_1.0_1)
surv_Sm_Ps_1.0_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_1.0_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30 30 18.3 7.51 19.9
## primary=S.m 25 25 36.7 3.74 19.9
##
## Chisq= 19.9 on 1 degrees of freedom, p= 8e-06
plot(survfit(Surv(Sm_Ps_block1$hour, Sm_Ps_block1$status) ~ Sm_Ps_block1$primary + Sm_Ps_block1$secondary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1), lwd=2.3, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Hours Post-Secondary Infection", cex.lab=1.3, xlim=c(10,48), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-PBS", "0.001", "PBS-0.01", "PBS-0.1", "PBS-1.0", "S.m-PBS", "S.m-0.001", "S.m-0.01", "S.m-0.1", "S.m-1.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1),lwd=4, cex=1)
Sm_Ps_block2 <- subset(Sm_Ps, date == "10/10/2020")
table(Sm_Ps_block2$primary, Sm_Ps_block2$secondary_dose)
##
## 0 0.001 0.01 0.1 1
## PBS 20 30 30 30 30
## S.m 20 30 30 30 30
#subset by secondary dose
Sm_Ps_0.001_2<-subset(Sm_Ps_block2, secondary_dose=="0.001")
Sm_Ps_0.01_2<-subset(Sm_Ps_block2, secondary_dose=="0.01")
Sm_Ps_0.1_2<-subset(Sm_Ps_block2, secondary_dose=="0.1")
Sm_Ps_1.0_2<-subset(Sm_Ps_block2, secondary_dose=="1")
surv_Sm_Ps_0.001_2<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_0.001_2)
surv_Sm_Ps_0.001_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_0.001_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30 25 16 5.07 10.5
## primary=S.m 30 13 22 3.69 10.5
##
## Chisq= 10.5 on 1 degrees of freedom, p= 0.001
surv_Sm_Ps_0.01_2<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_0.01_2)
surv_Sm_Ps_0.01_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_0.01_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30 30 15.8 12.77 24
## primary=S.m 30 26 40.2 5.02 24
##
## Chisq= 24 on 1 degrees of freedom, p= 9e-07
surv_Sm_Ps_0.1_2<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_0.1_2)
surv_Sm_Ps_0.1_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_0.1_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30 30 14.6 16.35 40.9
## primary=S.m 30 30 45.4 5.24 40.9
##
## Chisq= 40.9 on 1 degrees of freedom, p= 2e-10
surv_Sm_Ps_1.0_2<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_1.0_2)
surv_Sm_Ps_1.0_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_1.0_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30 30 13.5 20.14 44
## primary=S.m 30 30 46.5 5.85 44
##
## Chisq= 44 on 1 degrees of freedom, p= 3e-11
plot(survfit(Surv(Sm_Ps_block2$hour, Sm_Ps_block2$status) ~ Sm_Ps_block2$primary + Sm_Ps_block2$secondary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1), lwd=2.3, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Hours Post-Secondary Infection", cex.lab=1.3, xlim=c(10,48), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-PBS", "0.001", "PBS-0.01", "PBS-0.1", "PBS-1.0", "S.m-PBS", "S.m-0.001", "S.m-0.01", "S.m-0.1", "S.m-1.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1),lwd=4, cex=1)
Sm_Ps_block3 <- subset(Sm_Ps, date == "10/24/2020")
table(Sm_Ps_block3$primary, Sm_Ps_block3$secondary_dose)
##
## 0 0.001 0.01 0.1 1
## PBS 20 22 28 30 15
## S.m 19 20 26 22 26
#subset by secondary dose
Sm_Ps_0.001_3<-subset(Sm_Ps_block3, secondary_dose=="0.001")
Sm_Ps_0.01_3<-subset(Sm_Ps_block3, secondary_dose=="0.01")
Sm_Ps_0.1_3<-subset(Sm_Ps_block3, secondary_dose=="0.1")
Sm_Ps_1.0_3<-subset(Sm_Ps_block3, secondary_dose=="1")
surv_Sm_Ps_0.001_3<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_0.001_3)
surv_Sm_Ps_0.001_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_0.001_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 22 21 15.6 1.90 9.75
## primary=S.m 20 13 18.4 1.61 9.75
##
## Chisq= 9.8 on 1 degrees of freedom, p= 0.002
surv_Sm_Ps_0.01_3<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_0.01_3)
surv_Sm_Ps_0.01_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_0.01_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 28 28 13.1 16.80 34.7
## primary=S.m 26 21 35.9 6.16 34.7
##
## Chisq= 34.7 on 1 degrees of freedom, p= 4e-09
surv_Sm_Ps_0.1_3<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_0.1_3)
surv_Sm_Ps_0.1_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_0.1_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30 30 14 18.33 40.2
## primary=S.m 22 21 37 6.93 40.2
##
## Chisq= 40.2 on 1 degrees of freedom, p= 2e-10
surv_Sm_Ps_1.0_3<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_1.0_3)
surv_Sm_Ps_1.0_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_1.0_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 15 15 4.8 21.71 34.2
## primary=S.m 26 26 36.2 2.88 34.2
##
## Chisq= 34.2 on 1 degrees of freedom, p= 5e-09
plot(survfit(Surv(Sm_Ps_block3$hour, Sm_Ps_block3$status) ~ Sm_Ps_block3$primary + Sm_Ps_block3$secondary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1), lwd=2.3, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Hours Post-Secondary Infection", cex.lab=1.3, xlim=c(10,48), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-PBS", "0.001", "PBS-0.01", "PBS-0.1", "PBS-1.0", "S.m-PBS", "S.m-0.001", "S.m-0.01", "S.m-0.1", "S.m-1.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1),lwd=4, cex=1)
# Grab p-values from each individual logrank test (results shown above in each block)
pvalue_Sm_Ps_0.001_1<-pchisq(surv_Sm_Ps_0.001_1$chisq, df=1, lower.tail=F)
pvalue_Sm_Ps_0.01_1<-pchisq(surv_Sm_Ps_0.01_1$chisq, df=1, lower.tail=F)
pvalue_Sm_Ps_0.1_1<-pchisq(surv_Sm_Ps_0.1_1$chisq, df=1, lower.tail=F)
pvalue_Sm_Ps_1.0_1<-pchisq(surv_Sm_Ps_1.0_1$chisq, df=1, lower.tail=F)
pvalue_Sm_Ps_0.001_2<-pchisq(surv_Sm_Ps_0.001_2$chisq, df=1, lower.tail=F)
pvalue_Sm_Ps_0.01_2<-pchisq(surv_Sm_Ps_0.01_2$chisq, df=1, lower.tail=F)
pvalue_Sm_Ps_0.1_2<-pchisq(surv_Sm_Ps_0.1_2$chisq, df=1, lower.tail=F)
pvalue_Sm_Ps_1.0_2<-pchisq(surv_Sm_Ps_1.0_2$chisq, df=1, lower.tail=F)
pvalue_Sm_Ps_0.001_3<-pchisq(surv_Sm_Ps_0.001_3$chisq, df=1, lower.tail=F)
pvalue_Sm_Ps_0.01_3<-pchisq(surv_Sm_Ps_0.01_3$chisq, df=1, lower.tail=F)
pvalue_Sm_Ps_0.1_3<-pchisq(surv_Sm_Ps_0.1_3$chisq, df=1, lower.tail=F)
pvalue_Sm_Ps_1.0_3<-pchisq(surv_Sm_Ps_1.0_3$chisq, df=1, lower.tail=F)
# Use all p-values from a given condition to get Fisher's combined test p-value
pvalues_Sm_Ps_0.001 <- c(pvalue_Sm_Ps_0.001_2,pvalue_Sm_Ps_0.001_3)
test.stat_Sm_Ps_0.001 <- -2*sum(log(pvalues_Sm_Ps_0.001))
pfinal_Sm_Ps_0.001<-pchisq(test.stat_Sm_Ps_0.001, df=2*length(pvalues_Sm_Ps_0.001), lower.tail=F)
pvalues_Sm_Ps_0.01 <- c(pvalue_Sm_Ps_0.01_1,pvalue_Sm_Ps_0.01_2,pvalue_Sm_Ps_0.01_3)
test.stat_Sm_Ps_0.01 <- -2*sum(log(pvalues_Sm_Ps_0.01))
pfinal_Sm_Ps_0.01<-pchisq(test.stat_Sm_Ps_0.01, df=2*length(pvalues_Sm_Ps_0.01), lower.tail=F)
pvalues_Sm_Ps_0.1 <- c(pvalue_Sm_Ps_0.1_1,pvalue_Sm_Ps_0.1_2,pvalue_Sm_Ps_0.1_3)
test.stat_Sm_Ps_0.1 <- -2*sum(log(pvalues_Sm_Ps_0.1))
pfinal_Sm_Ps_0.1<-pchisq(test.stat_Sm_Ps_0.1, df=2*length(pvalues_Sm_Ps_0.1), lower.tail=F)
pvalues_Sm_Ps_1.0 <- c(pvalue_Sm_Ps_1.0_1,pvalue_Sm_Ps_1.0_2,pvalue_Sm_Ps_1.0_3)
test.stat_Sm_Ps_1.0 <- -2*sum(log(pvalues_Sm_Ps_1.0))
pfinal_Sm_Ps_1.0<-pchisq(test.stat_Sm_Ps_1.0, df=2*length(pvalues_Sm_Ps_1.0), lower.tail=F)
pvalues_Sm_Ps_0.001_table <- c(pvalue_Sm_Ps_0.001_1, pvalue_Sm_Ps_0.001_2, pvalue_Sm_Ps_0.001_3)
p_Sm_Ps_0.001<-scales::pvalue(pvalues_Sm_Ps_0.001_table)
pvalues_Sm_Ps_0.01_table <- c(pvalue_Sm_Ps_0.01_1,pvalue_Sm_Ps_0.01_2,pvalue_Sm_Ps_0.01_3)
p_Sm_Ps_0.01<-scales::pvalue(pvalues_Sm_Ps_0.01_table)
pvalues_Sm_Ps_0.1_table <- c(pvalue_Sm_Ps_0.1_1,pvalue_Sm_Ps_0.1_2,pvalue_Sm_Ps_0.1_3)
p_Sm_Ps_0.1<-scales::pvalue(pvalues_Sm_Ps_0.1_table)
pvalues_Sm_Ps_1.0_table <- c(pvalue_Sm_Ps_1.0_1,pvalue_Sm_Ps_1.0_2,pvalue_Sm_Ps_1.0_3)
p_Sm_Ps_1.0<-scales::pvalue(pvalues_Sm_Ps_1.0_table)
pvalue_table_Sm_Ps <- data.frame("0.001" = p_Sm_Ps_0.001, "0.01" = p_Sm_Ps_0.01, "0.1" = p_Sm_Ps_0.1, "1.0" = p_Sm_Ps_1.0)
colnames(pvalue_table_Sm_Ps) <- c(0.001,0.01, 0.1, 1.0)
pfisher_Sm_Ps<- c(pfinal_Sm_Ps_0.001, pfinal_Sm_Ps_0.01, pfinal_Sm_Ps_0.1, pfinal_Sm_Ps_1.0)
pfisher_Sm_Ps<-scales::pvalue(pfisher_Sm_Ps)
pvalue_table_Sm_Ps[nrow(pvalue_table_Sm_Ps) + 1,]<-pfisher_Sm_Ps
row.names(pvalue_table_Sm_Ps) <- c("Block 1 - October 1st, 2021", "Block 2 - October 10th, 2021", "Block 3 - October 24th, 2021", " Final P-value ")
pvalue_table_Sm_Ps %>%
kbl(caption = "<center><strong>P-values for log-rank tests on individual doses within each experiment<center><strong>") %>%
kable_classic(full_width = F, html_font = "Cambria", position = "center")%>%
add_header_above(c(" "=1, "Infectious Dose (Abs600nm)" = 4)) %>%
pack_rows("Individual Blocks", 1, 3) %>%
pack_rows("Fisher's Combined", 4, 4)
|
Infectious Dose (Abs600nm)
|
||||
|---|---|---|---|---|
| 0.001 | 0.01 | 0.1 | 1 | |
| Individual Blocks | ||||
| Block 1 - October 1st, 2021 | 0.287 | <0.001 | <0.001 | <0.001 |
| Block 2 - October 10th, 2021 | 0.001 | <0.001 | <0.001 | <0.001 |
| Block 3 - October 24th, 2021 | 0.002 | <0.001 | <0.001 | <0.001 |
| Fisher’s Combined | ||||
| Final P-value | <0.001 | <0.001 | <0.001 | <0.001 |
Study design: + Initially
analyzed using Cox Proportional Hazards mixed effects model with date of
injection, primary injection dose as fixed effects and random effect
terms accounting for experimental structure (e.g. date of injection,
housing of groups of flies in vials). However, this model violated the
assumptions of Cox Proportional Hazards. + Therefore we switched to
using log-rank pairwise comparisons to determine whether dose used to
induce chronic infection significantly impacted survival of secondary
infection by comparing to the control flies given only a secondary
infection. + LogRank comparisons were done on data from each date
individually + If experiments on all dates showed higher survival in
conditions with chronic infection, p-values were combined using Fisher’s
test and only the final p-value is reported in the text. + If there was
a variable effect of chronic infection, the number of individual
experiments with significant effects is reported. +
Conclusions + Moderate doses of S. marcescens
(Abs600nm = 0.01, 0.1 & 1.0) all provided significant protection
against secondary infection with both P. rettgeri ad P.
sneebia (P<0.001) + The lowest dose of S. marcescens
(Abs600nm = 0.01, 0.1 & 1.0) did not provide protection against
either secondary infection. + The highest dose provided inconsistent
protection, sometimes providing robust protection and other times dying
as quickly as (or even faster than) control flies carrying only a
secondary infection.
Pr_varied_primary = subset(Pr_varied_primary, secondary == "P.rett")
Pr_varied_primary$primary_dose <- as.factor(Pr_varied_primary$primary_dose)
table(Pr_varied_primary$date, Pr_varied_primary$primary_dose)
##
## 0 0.001 0.01 0.1 1 2
## 10/18/2020 40 38 35 40 27 31
## 10/20/2020 36 37 35 39 35 30
## 10/29/2020 40 37 46 36 37 33
## 11/2/2020 37 31 30 40 30 27
## 11/28/2020 39 47 47 48 45 42
## 11/7/2020 40 38 40 35 34 36
## 12/4/2020 35 38 36 34 27 0
## 6/17/2021 37 47 47 44 42 36
## 6/18/2021 33 48 48 42 42 63
table(Pr_varied_primary$primary_dose)
##
## 0 0.001 0.01 0.1 1 2
## 337 361 364 358 319 298
#png(file="Figure5A.png", width = 1200, height = 900, res = 130)
par(mar=c(6,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(Pr_varied_primary$day, Pr_varied_primary$status) ~ Pr_varied_primary$primary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=5, yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i', ylim=c(0, 1), yaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=4.5)
box(lwd=5)
#legend("bottomleft", bty="n" , c("secondary infection only", "chronic infection (0.001)", "chronic infection (0.01)", "chronic infection (0.1)", "chronic infection (1.0)", "chronic infection (2.0)"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1.5)
#dev.off()
Pr_varied_primary_block1 <- subset(Pr_varied_primary, date == "10/18/2020")
table(Pr_varied_primary_block1$primary_dose)
##
## 0 0.001 0.01 0.1 1 2
## 40 38 35 40 27 31
#subset by primary dose for pairwise comparison dose
Pr_varied_primary_0.001_1<-subset(Pr_varied_primary_block1, primary_dose == "0.001"|primary_dose=="0")
table(Pr_varied_primary_0.001_1$primary_dose)
##
## 0 0.001 0.01 0.1 1 2
## 40 38 0 0 0 0
Pr_varied_primary_0.01_1<-subset(Pr_varied_primary_block1, primary_dose=="0.01"|primary_dose=="0")
Pr_varied_primary_0.1_1<-subset(Pr_varied_primary_block1, primary_dose=="0.1"|primary_dose=="0")
Pr_varied_primary_1.0_1<-subset(Pr_varied_primary_block1, primary_dose=="1"|primary_dose=="0")
Pr_varied_primary_2.0_1<-subset(Pr_varied_primary_block1, primary_dose=="2"|primary_dose=="0")
surv_Pr_varied_primary_0.001_1<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.001_1)
surv_Pr_varied_primary_0.001_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.001_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40 32 37.4 0.792 3.24
## primary=S.m 38 35 29.6 1.003 3.24
##
## Chisq= 3.2 on 1 degrees of freedom, p= 0.07
surv_Pr_varied_primary_0.01_1<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.01_1)
surv_Pr_varied_primary_0.01_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.01_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40 32 28.8 0.358 1.11
## primary=S.m 35 24 27.2 0.379 1.11
##
## Chisq= 1.1 on 1 degrees of freedom, p= 0.3
surv_Pr_varied_primary_0.1_1<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.1_1)
surv_Pr_varied_primary_0.1_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.1_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40 32 18.6 9.60 22.5
## primary=S.m 40 12 25.4 7.05 22.5
##
## Chisq= 22.5 on 1 degrees of freedom, p= 2e-06
surv_Pr_varied_primary_1.0_1<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_1.0_1)
surv_Pr_varied_primary_1.0_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_1.0_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40 32 26.3 1.22 3.27
## primary=S.m 27 22 27.7 1.16 3.27
##
## Chisq= 3.3 on 1 degrees of freedom, p= 0.07
surv_Pr_varied_primary_2.0_1<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_2.0_1)
surv_Pr_varied_primary_2.0_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_2.0_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40 32 37.4 0.776 4.16
## primary=S.m 31 31 25.6 1.133 4.16
##
## Chisq= 4.2 on 1 degrees of freedom, p= 0.04
par(mar=c(4,4,0.1,0.1))
plot(survfit(Surv(Pr_varied_primary_block1$day, Pr_varied_primary_block1$status) ~ Pr_varied_primary_block1$primary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Days Post-Secondary Infection", cex.lab=1.3, xlim=c(0,7), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-P.rett 2.0", "S.m 0.001-P.rett 2.0", "S.m 0.01-P.rett 2.0", "S.m 0.1-P.rett 2.0", "S.m 1.0-P.rett 2.0", "S.m 2.0-P.rett 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)
Pr_varied_primary_block2 <- subset(Pr_varied_primary, date == "10/20/2020")
table(Pr_varied_primary_block2$primary_dose)
##
## 0 0.001 0.01 0.1 1 2
## 36 37 35 39 35 30
#subset by primary dose for pairwise comparison dose
Pr_varied_primary_0.001_2<-subset(Pr_varied_primary_block2, primary_dose == "0.001"|primary_dose=="0")
Pr_varied_primary_0.01_2<-subset(Pr_varied_primary_block2, primary_dose=="0.01"|primary_dose=="0")
Pr_varied_primary_0.1_2<-subset(Pr_varied_primary_block2, primary_dose=="0.1"|primary_dose=="0")
Pr_varied_primary_1.0_2<-subset(Pr_varied_primary_block2, primary_dose=="1"|primary_dose=="0")
Pr_varied_primary_2.0_2<-subset(Pr_varied_primary_block2, primary_dose=="2"|primary_dose=="0")
surv_Pr_varied_primary_0.001_2<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.001_2)
surv_Pr_varied_primary_0.001_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.001_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 36 31 32.6 0.0783 0.403
## primary=S.m 37 35 33.4 0.0764 0.403
##
## Chisq= 0.4 on 1 degrees of freedom, p= 0.5
surv_Pr_varied_primary_0.01_2<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.01_2)
surv_Pr_varied_primary_0.01_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.01_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 36 31 27.7 0.388 1.57
## primary=S.m 35 28 31.3 0.344 1.57
##
## Chisq= 1.6 on 1 degrees of freedom, p= 0.2
surv_Pr_varied_primary_0.1_2<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.1_2)
surv_Pr_varied_primary_0.1_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.1_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 36 31 19.6 6.63 16.7
## primary=S.m 39 18 29.4 4.42 16.7
##
## Chisq= 16.7 on 1 degrees of freedom, p= 4e-05
surv_Pr_varied_primary_1.0_2<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_1.0_2)
surv_Pr_varied_primary_1.0_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_1.0_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 36 31 21.7 4.03 9.17
## primary=S.m 35 28 37.3 2.34 9.17
##
## Chisq= 9.2 on 1 degrees of freedom, p= 0.002
surv_Pr_varied_primary_2.0_2<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_2.0_2)
surv_Pr_varied_primary_2.0_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_2.0_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 36 31 25.5 1.192 3.23
## primary=S.m 30 27 32.5 0.935 3.23
##
## Chisq= 3.2 on 1 degrees of freedom, p= 0.07
par(mar=c(4,4,0.1,0.1))
plot(survfit(Surv(Pr_varied_primary_block2$day, Pr_varied_primary_block2$status) ~ Pr_varied_primary_block2$primary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Days Post-Secondary Infection", cex.lab=1.3, xlim=c(0,7), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-P.rett 2.0", "S.m 0.001-P.rett 2.0", "S.m 0.01-P.rett 2.0", "S.m 0.1-P.rett 2.0", "S.m 1.0-P.rett 2.0", "S.m 2.0-P.rett 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)
Pr_varied_primary_block3 <- subset(Pr_varied_primary, date == "10/29/2020")
table(Pr_varied_primary_block3$primary_dose)
##
## 0 0.001 0.01 0.1 1 2
## 40 37 46 36 37 33
#subset by primary dose for pairwise comparison dose
Pr_varied_primary_0.001_3<-subset(Pr_varied_primary_block3, primary_dose == "0.001"|primary_dose=="0")
Pr_varied_primary_0.01_3<-subset(Pr_varied_primary_block3, primary_dose=="0.01"|primary_dose=="0")
Pr_varied_primary_0.1_3<-subset(Pr_varied_primary_block3, primary_dose=="0.1"|primary_dose=="0")
Pr_varied_primary_1.0_3<-subset(Pr_varied_primary_block3, primary_dose=="1"|primary_dose=="0")
Pr_varied_primary_2.0_3<-subset(Pr_varied_primary_block3, primary_dose=="2"|primary_dose=="0")
surv_Pr_varied_primary_0.001_3<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.001_3)
surv_Pr_varied_primary_0.001_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.001_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40 36 34.6 0.0559 0.235
## primary=S.m 37 31 32.4 0.0597 0.235
##
## Chisq= 0.2 on 1 degrees of freedom, p= 0.6
surv_Pr_varied_primary_0.01_3<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.01_3)
surv_Pr_varied_primary_0.01_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.01_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40 36 28.6 1.93 5.61
## primary=S.m 46 32 39.4 1.40 5.61
##
## Chisq= 5.6 on 1 degrees of freedom, p= 0.02
surv_Pr_varied_primary_0.1_3<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.1_3)
surv_Pr_varied_primary_0.1_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.1_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40 36 18.3 17.1 41.4
## primary=S.m 36 9 26.7 11.7 41.4
##
## Chisq= 41.4 on 1 degrees of freedom, p= 1e-10
surv_Pr_varied_primary_1.0_3<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_1.0_3)
surv_Pr_varied_primary_1.0_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_1.0_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40 36 19.6 13.83 30.3
## primary=S.m 37 19 35.4 7.63 30.3
##
## Chisq= 30.3 on 1 degrees of freedom, p= 4e-08
surv_Pr_varied_primary_2.0_3<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_2.0_3)
surv_Pr_varied_primary_2.0_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_2.0_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40 36 27.8 2.43 6.33
## primary=S.m 33 27 35.2 1.92 6.33
##
## Chisq= 6.3 on 1 degrees of freedom, p= 0.01
par(mar=c(4,4,0.1,0.1))
plot(survfit(Surv(Pr_varied_primary_block3$day, Pr_varied_primary_block3$status) ~ Pr_varied_primary_block3$primary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Days Post-Secondary Infection", cex.lab=1.3, xlim=c(0,7), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-P.rett 2.0", "S.m 0.001-P.rett 2.0", "S.m 0.01-P.rett 2.0", "S.m 0.1-P.rett 2.0", "S.m 1.0-P.rett 2.0", "S.m 2.0-P.rett 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)
Pr_varied_primary_block4 <- subset(Pr_varied_primary, date == "11/2/2020")
table(Pr_varied_primary_block4$primary_dose)
##
## 0 0.001 0.01 0.1 1 2
## 37 31 30 40 30 27
#subset by primary dose for pairwise comparison dose
Pr_varied_primary_0.001_4<-subset(Pr_varied_primary_block4, primary_dose == "0.001"|primary_dose=="0")
Pr_varied_primary_0.01_4<-subset(Pr_varied_primary_block4, primary_dose=="0.01"|primary_dose=="0")
Pr_varied_primary_0.1_4<-subset(Pr_varied_primary_block4, primary_dose=="0.1"|primary_dose=="0")
Pr_varied_primary_1.0_4<-subset(Pr_varied_primary_block4, primary_dose=="1"|primary_dose=="0")
Pr_varied_primary_2.0_4<-subset(Pr_varied_primary_block4, primary_dose=="2"|primary_dose=="0")
surv_Pr_varied_primary_0.001_4<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.001_4)
surv_Pr_varied_primary_0.001_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.001_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 33 33.2 0.000817 0.00537
## primary=S.m 31 28 27.8 0.000973 0.00537
##
## Chisq= 0 on 1 degrees of freedom, p= 0.9
surv_Pr_varied_primary_0.01_4<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.01_4)
surv_Pr_varied_primary_0.01_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.01_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 33 26.5 1.57 6.91
## primary=S.m 30 20 26.5 1.58 6.91
##
## Chisq= 6.9 on 1 degrees of freedom, p= 0.009
surv_Pr_varied_primary_0.1_4<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.1_4)
surv_Pr_varied_primary_0.1_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.1_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 33 21.2 6.58 20.5
## primary=S.m 40 19 30.8 4.52 20.5
##
## Chisq= 20.5 on 1 degrees of freedom, p= 6e-06
surv_Pr_varied_primary_1.0_4<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_1.0_4)
surv_Pr_varied_primary_1.0_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_1.0_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 33 20.9 7.08 19.2
## primary=S.m 30 19 31.1 4.74 19.2
##
## Chisq= 19.2 on 1 degrees of freedom, p= 1e-05
surv_Pr_varied_primary_2.0_4<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_2.0_4)
surv_Pr_varied_primary_2.0_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_2.0_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 33 30.6 0.181 0.686
## primary=S.m 27 27 29.4 0.189 0.686
##
## Chisq= 0.7 on 1 degrees of freedom, p= 0.4
par(mar=c(4,4,0.1,0.1))
plot(survfit(Surv(Pr_varied_primary_block4$day, Pr_varied_primary_block4$status) ~ Pr_varied_primary_block4$primary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Days Post-Secondary Infection", cex.lab=1.3, xlim=c(0,7), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-P.rett 2.0", "S.m 0.001-P.rett 2.0", "S.m 0.01-P.rett 2.0", "S.m 0.1-P.rett 2.0", "S.m 1.0-P.rett 2.0", "S.m 2.0-P.rett 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)
Pr_varied_primary_block5 <- subset(Pr_varied_primary, date == "11/7/2020")
table(Pr_varied_primary_block5$primary_dose)
##
## 0 0.001 0.01 0.1 1 2
## 40 38 40 35 34 36
#subset by primary dose for pairwise comparison dose
Pr_varied_primary_0.001_5<-subset(Pr_varied_primary_block5, primary_dose == "0.001"|primary_dose=="0")
Pr_varied_primary_0.01_5<-subset(Pr_varied_primary_block5, primary_dose=="0.01"|primary_dose=="0")
Pr_varied_primary_0.1_5<-subset(Pr_varied_primary_block5, primary_dose=="0.1"|primary_dose=="0")
Pr_varied_primary_1.0_5<-subset(Pr_varied_primary_block5, primary_dose=="1"|primary_dose=="0")
Pr_varied_primary_2.0_5<-subset(Pr_varied_primary_block5, primary_dose=="2"|primary_dose=="0")
surv_Pr_varied_primary_0.001_5<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.001_5)
surv_Pr_varied_primary_0.001_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.001_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40 35 33.9 0.0383 0.156
## primary=S.m 38 32 33.1 0.0391 0.156
##
## Chisq= 0.2 on 1 degrees of freedom, p= 0.7
surv_Pr_varied_primary_0.01_5<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.01_5)
surv_Pr_varied_primary_0.01_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.01_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40 35 32 0.289 0.942
## primary=S.m 40 35 38 0.243 0.942
##
## Chisq= 0.9 on 1 degrees of freedom, p= 0.3
surv_Pr_varied_primary_0.1_5<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.1_5)
surv_Pr_varied_primary_0.1_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.1_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40 35 21.5 8.41 21.9
## primary=S.m 35 15 28.5 6.37 21.9
##
## Chisq= 21.9 on 1 degrees of freedom, p= 3e-06
surv_Pr_varied_primary_1.0_5<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_1.0_5)
surv_Pr_varied_primary_1.0_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_1.0_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40 35 20 11.28 27.4
## primary=S.m 34 16 31 7.27 27.4
##
## Chisq= 27.4 on 1 degrees of freedom, p= 2e-07
surv_Pr_varied_primary_2.0_5<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_2.0_5)
surv_Pr_varied_primary_2.0_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_2.0_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40 35 26.1 3.02 6.7
## primary=S.m 36 34 42.9 1.84 6.7
##
## Chisq= 6.7 on 1 degrees of freedom, p= 0.01
par(mar=c(4,4,0.1,0.1))
plot(survfit(Surv(Pr_varied_primary_block5$day, Pr_varied_primary_block5$status) ~ Pr_varied_primary_block5$primary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Days Post-Secondary Infection", cex.lab=1.3, xlim=c(0,7), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-P.rett 2.0", "S.m 0.001-P.rett 2.0", "S.m 0.01-P.rett 2.0", "S.m 0.1-P.rett 2.0", "S.m 1.0-P.rett 2.0", "S.m 2.0-P.rett 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)
Pr_varied_primary_block6 <- subset(Pr_varied_primary, date == "11/28/2020")
table(Pr_varied_primary_block6$primary_dose)
##
## 0 0.001 0.01 0.1 1 2
## 39 47 47 48 45 42
#subset by primary dose for pairwise comparison dose
Pr_varied_primary_0.001_6<-subset(Pr_varied_primary_block6, primary_dose == "0.001"|primary_dose=="0")
Pr_varied_primary_0.01_6<-subset(Pr_varied_primary_block6, primary_dose=="0.01"|primary_dose=="0")
Pr_varied_primary_0.1_6<-subset(Pr_varied_primary_block6, primary_dose=="0.1"|primary_dose=="0")
Pr_varied_primary_1.0_6<-subset(Pr_varied_primary_block6, primary_dose=="1"|primary_dose=="0")
Pr_varied_primary_2.0_6<-subset(Pr_varied_primary_block6, primary_dose=="2"|primary_dose=="0")
surv_Pr_varied_primary_0.001_6<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.001_6)
surv_Pr_varied_primary_0.001_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.001_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 39 33 38.1 0.679 3.44
## primary=S.m 47 44 38.9 0.664 3.44
##
## Chisq= 3.4 on 1 degrees of freedom, p= 0.06
surv_Pr_varied_primary_0.01_6<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.01_6)
surv_Pr_varied_primary_0.01_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.01_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 39 33 32.7 0.00346 0.0119
## primary=S.m 47 39 39.3 0.00288 0.0119
##
## Chisq= 0 on 1 degrees of freedom, p= 0.9
surv_Pr_varied_primary_0.1_6<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.1_6)
surv_Pr_varied_primary_0.1_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.1_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 39 33 19.6 9.14 19.7
## primary=S.m 48 23 36.4 4.93 19.7
##
## Chisq= 19.7 on 1 degrees of freedom, p= 9e-06
surv_Pr_varied_primary_1.0_6<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_1.0_6)
surv_Pr_varied_primary_1.0_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_1.0_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 39 33 19.1 10.03 21.9
## primary=S.m 45 20 33.9 5.67 21.9
##
## Chisq= 21.9 on 1 degrees of freedom, p= 3e-06
surv_Pr_varied_primary_2.0_6<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_2.0_6)
surv_Pr_varied_primary_2.0_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_2.0_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 39 33 25.8 2.03 4.56
## primary=S.m 42 37 44.2 1.18 4.56
##
## Chisq= 4.6 on 1 degrees of freedom, p= 0.03
par(mar=c(4,4,0.1,0.1))
plot(survfit(Surv(Pr_varied_primary_block6$day, Pr_varied_primary_block6$status) ~ Pr_varied_primary_block6$primary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Days Post-Secondary Infection", cex.lab=1.3, xlim=c(0,7), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-P.rett 2.0", "S.m 0.001-P.rett 2.0", "S.m 0.01-P.rett 2.0", "S.m 0.1-P.rett 2.0", "S.m 1.0-P.rett 2.0", "S.m 2.0-P.rett 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)
Pr_varied_primary_block7 <- subset(Pr_varied_primary, date == "12/4/2020")
table(Pr_varied_primary_block7$primary_dose)
##
## 0 0.001 0.01 0.1 1 2
## 35 38 36 34 27 0
#subset by primary dose for pairwise comparison dose
Pr_varied_primary_0.001_7<-subset(Pr_varied_primary_block7, primary_dose == "0.001"|primary_dose=="0")
Pr_varied_primary_0.01_7<-subset(Pr_varied_primary_block7, primary_dose=="0.01"|primary_dose=="0")
Pr_varied_primary_0.1_7<-subset(Pr_varied_primary_block7, primary_dose=="0.1"|primary_dose=="0")
Pr_varied_primary_1.0_7<-subset(Pr_varied_primary_block7, primary_dose=="1"|primary_dose=="0")
#Pr_varied_primary_2.0_7<-subset(Pr_varied_primary_block7, primary_dose=="2"|primary_dose=="0")
surv_Pr_varied_primary_0.001_7<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.001_7)
surv_Pr_varied_primary_0.001_7
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.001_7)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 35 31 31.1 0.000230 0.00162
## primary=S.m 38 33 32.9 0.000217 0.00162
##
## Chisq= 0 on 1 degrees of freedom, p= 1
surv_Pr_varied_primary_0.01_7<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.01_7)
surv_Pr_varied_primary_0.01_7
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.01_7)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 35 31 24.8 1.53 7.31
## primary=S.m 36 21 27.2 1.40 7.31
##
## Chisq= 7.3 on 1 degrees of freedom, p= 0.007
surv_Pr_varied_primary_0.1_7<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.1_7)
surv_Pr_varied_primary_0.1_7
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.1_7)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 35 31 16.1 13.8 40.8
## primary=S.m 34 6 20.9 10.6 40.8
##
## Chisq= 40.8 on 1 degrees of freedom, p= 2e-10
surv_Pr_varied_primary_1.0_7<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_1.0_7)
surv_Pr_varied_primary_1.0_7
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_1.0_7)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 35 31 20.8 4.96 14.1
## primary=S.m 27 21 31.2 3.31 14.1
##
## Chisq= 14.1 on 1 degrees of freedom, p= 2e-04
#surv_Pr_varied_primary_2.0_7<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_2.0_7)
#surv_Pr_varied_primary_2.0_7
par(mar=c(4,4,0.1,0.1))
plot(survfit(Surv(Pr_varied_primary_block7$day, Pr_varied_primary_block7$status) ~ Pr_varied_primary_block7$primary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Days Post-Secondary Infection", cex.lab=1.3, xlim=c(0,7), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-P.rett 2.0", "S.m 0.001-P.rett 2.0", "S.m 0.01-P.rett 2.0", "S.m 0.1-P.rett 2.0", "S.m 1.0-P.rett 2.0", "S.m 2.0-P.rett 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)
Pr_varied_primary_block8 <- subset(Pr_varied_primary, date == "6/17/2021")
table(Pr_varied_primary_block8$primary_dose)
##
## 0 0.001 0.01 0.1 1 2
## 37 47 47 44 42 36
#subset by primary dose for pairwise comparison dose
Pr_varied_primary_0.001_8<-subset(Pr_varied_primary_block8, primary_dose == "0.001"|primary_dose=="0")
Pr_varied_primary_0.01_8<-subset(Pr_varied_primary_block8, primary_dose=="0.01"|primary_dose=="0")
Pr_varied_primary_0.1_8<-subset(Pr_varied_primary_block8, primary_dose=="0.1"|primary_dose=="0")
Pr_varied_primary_1.0_8<-subset(Pr_varied_primary_block8, primary_dose=="1"|primary_dose=="0")
Pr_varied_primary_2.0_8<-subset(Pr_varied_primary_block8, primary_dose=="2"|primary_dose=="0")
surv_Pr_varied_primary_0.001_8<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.001_8)
surv_Pr_varied_primary_0.001_8
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.001_8)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 37 36.2 0.0185 0.123
## primary=S.m 47 44 44.8 0.0150 0.123
##
## Chisq= 0.1 on 1 degrees of freedom, p= 0.7
surv_Pr_varied_primary_0.01_8<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.01_8)
surv_Pr_varied_primary_0.01_8
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.01_8)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 37 32.2 0.700 3.45
## primary=S.m 47 42 46.8 0.483 3.45
##
## Chisq= 3.4 on 1 degrees of freedom, p= 0.06
surv_Pr_varied_primary_0.1_8<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.1_8)
surv_Pr_varied_primary_0.1_8
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.1_8)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 37 17.5 21.9 55.4
## primary=S.m 44 14 33.5 11.4 55.4
##
## Chisq= 55.4 on 1 degrees of freedom, p= 1e-13
surv_Pr_varied_primary_1.0_8<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_1.0_8)
surv_Pr_varied_primary_1.0_8
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_1.0_8)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 37 16 27.4 67.1
## primary=S.m 42 13 34 12.9 67.1
##
## Chisq= 67.1 on 1 degrees of freedom, p= 3e-16
surv_Pr_varied_primary_2.0_8<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_2.0_8)
surv_Pr_varied_primary_2.0_8
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_2.0_8)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 37 24 7.02 27.6
## primary=S.m 36 35 48 3.51 27.6
##
## Chisq= 27.6 on 1 degrees of freedom, p= 2e-07
par(mar=c(4,4,0.1,0.1))
plot(survfit(Surv(Pr_varied_primary_block8$day, Pr_varied_primary_block8$status) ~ Pr_varied_primary_block8$primary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Days Post-Secondary Infection", cex.lab=1.3, xlim=c(0,7), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-P.rett 2.0", "S.m 0.001-P.rett 2.0", "S.m 0.01-P.rett 2.0", "S.m 0.1-P.rett 2.0", "S.m 1.0-P.rett 2.0", "S.m 2.0-P.rett 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)
Pr_varied_primary_block9<- subset(Pr_varied_primary, date == "6/18/2021")
table(Pr_varied_primary_block9$primary_dose)
##
## 0 0.001 0.01 0.1 1 2
## 33 48 48 42 42 63
#subset by primary dose for pairwise comparison dose
Pr_varied_primary_0.001_9<-subset(Pr_varied_primary_block9, primary_dose == "0.001"|primary_dose=="0")
Pr_varied_primary_0.01_9<-subset(Pr_varied_primary_block9, primary_dose=="0.01"|primary_dose=="0")
Pr_varied_primary_0.1_9<-subset(Pr_varied_primary_block9, primary_dose=="0.1"|primary_dose=="0")
Pr_varied_primary_1.0_9<-subset(Pr_varied_primary_block9, primary_dose=="1"|primary_dose=="0")
Pr_varied_primary_2.0_9<-subset(Pr_varied_primary_block9, primary_dose=="2"|primary_dose=="0")
surv_Pr_varied_primary_0.001_9<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.001_9)
surv_Pr_varied_primary_0.001_9
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.001_9)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 33 33 32.2 0.0206 1.39
## primary=S.m 48 46 46.8 0.0142 1.39
##
## Chisq= 1.4 on 1 degrees of freedom, p= 0.2
surv_Pr_varied_primary_0.01_9<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.01_9)
surv_Pr_varied_primary_0.01_9
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.01_9)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 33 33 27.7 1.013 10.5
## primary=S.m 48 38 43.3 0.648 10.5
##
## Chisq= 10.5 on 1 degrees of freedom, p= 0.001
surv_Pr_varied_primary_0.1_9<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.1_9)
surv_Pr_varied_primary_0.1_9
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.1_9)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 33 33 19.8 8.80 38.8
## primary=S.m 42 18 31.2 5.58 38.8
##
## Chisq= 38.8 on 1 degrees of freedom, p= 5e-10
surv_Pr_varied_primary_1.0_9<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_1.0_9)
surv_Pr_varied_primary_1.0_9
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_1.0_9)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 33 33 17.2 14.6 53.7
## primary=S.m 42 13 28.8 8.7 53.7
##
## Chisq= 53.7 on 1 degrees of freedom, p= 2e-13
surv_Pr_varied_primary_2.0_9<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_2.0_9)
surv_Pr_varied_primary_2.0_9
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_2.0_9)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 33 33 19.2 9.82 35.5
## primary=S.m 63 50 63.8 2.97 35.5
##
## Chisq= 35.5 on 1 degrees of freedom, p= 2e-09
par(mar=c(4,4,0.1,0.1))
plot(survfit(Surv(Pr_varied_primary_block9$day, Pr_varied_primary_block9$status) ~ Pr_varied_primary_block9$primary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Days Post-Secondary Infection", cex.lab=1.3, xlim=c(0,7), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-P.rett 2.0", "S.m 0.001-P.rett 2.0", "S.m 0.01-P.rett 2.0", "S.m 0.1-P.rett 2.0", "S.m 1.0-P.rett 2.0", "S.m 2.0-P.rett 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)
# Grab p-values from each individual logrank test (results shown above in each block)
pvalue_Pr_varied_primary_0.001_1<-pchisq(surv_Pr_varied_primary_0.001_1$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.01_1<-pchisq(surv_Pr_varied_primary_0.01_1$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.1_1<-pchisq(surv_Pr_varied_primary_0.1_1$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_1.0_1<-pchisq(surv_Pr_varied_primary_1.0_1$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_2.0_1<-pchisq(surv_Pr_varied_primary_2.0_1$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.001_2<-pchisq(surv_Pr_varied_primary_0.001_2$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.01_2<-pchisq(surv_Pr_varied_primary_0.01_2$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.1_2<-pchisq(surv_Pr_varied_primary_0.1_2$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_1.0_2<-pchisq(surv_Pr_varied_primary_1.0_2$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_2.0_2<-pchisq(surv_Pr_varied_primary_2.0_2$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.001_3<-pchisq(surv_Pr_varied_primary_0.001_3$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.01_3<-pchisq(surv_Pr_varied_primary_0.01_3$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.1_3<-pchisq(surv_Pr_varied_primary_0.1_3$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_1.0_3<-pchisq(surv_Pr_varied_primary_1.0_3$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_2.0_3<-pchisq(surv_Pr_varied_primary_2.0_3$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.001_4<-pchisq(surv_Pr_varied_primary_0.001_4$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.01_4<-pchisq(surv_Pr_varied_primary_0.01_4$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.1_4<-pchisq(surv_Pr_varied_primary_0.1_4$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_1.0_4<-pchisq(surv_Pr_varied_primary_1.0_4$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_2.0_4<-pchisq(surv_Pr_varied_primary_2.0_4$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.001_5<-pchisq(surv_Pr_varied_primary_0.001_5$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.01_5<-pchisq(surv_Pr_varied_primary_0.01_5$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.1_5<-pchisq(surv_Pr_varied_primary_0.1_5$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_1.0_5<-pchisq(surv_Pr_varied_primary_1.0_5$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_2.0_5<-pchisq(surv_Pr_varied_primary_2.0_5$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.001_6<-pchisq(surv_Pr_varied_primary_0.001_6$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.01_6<-pchisq(surv_Pr_varied_primary_0.01_6$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.1_6<-pchisq(surv_Pr_varied_primary_0.1_6$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_1.0_6<-pchisq(surv_Pr_varied_primary_1.0_6$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_2.0_6<-pchisq(surv_Pr_varied_primary_2.0_6$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.001_7<-pchisq(surv_Pr_varied_primary_0.001_7$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.01_7<-pchisq(surv_Pr_varied_primary_0.01_7$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.1_7<-pchisq(surv_Pr_varied_primary_0.1_7$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_1.0_7<-pchisq(surv_Pr_varied_primary_1.0_7$chisq, df=1, lower.tail=F)
#pvalue_Pr_varied_primary_2.0_7<-pchisq(surv_Pr_varied_primary_2.0_7$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.001_8<-pchisq(surv_Pr_varied_primary_0.001_8$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.01_8<-pchisq(surv_Pr_varied_primary_0.01_8$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.1_8<-pchisq(surv_Pr_varied_primary_0.1_8$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_1.0_8<-pchisq(surv_Pr_varied_primary_1.0_8$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_2.0_8<-pchisq(surv_Pr_varied_primary_2.0_8$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.001_9<-pchisq(surv_Pr_varied_primary_0.001_9$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.01_9<-pchisq(surv_Pr_varied_primary_0.01_9$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.1_9<-pchisq(surv_Pr_varied_primary_0.1_9$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_1.0_9<-pchisq(surv_Pr_varied_primary_1.0_9$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_2.0_9<-pchisq(surv_Pr_varied_primary_2.0_9$chisq, df=1, lower.tail=F)
# Use all p-values from a given condition to get Fisher's combined test p-value
pvalues_Pr_varied_primary_0.001 <- c(pvalue_Pr_varied_primary_0.001_2 ,pvalue_Pr_varied_primary_0.001_3, pvalue_Pr_varied_primary_0.001_4, pvalue_Pr_varied_primary_0.001_5, pvalue_Pr_varied_primary_0.001_6, pvalue_Pr_varied_primary_0.001_7, pvalue_Pr_varied_primary_0.001_8, pvalue_Pr_varied_primary_0.001_9)
test.stat_Pr_varied_primary_0.001 <- -2*sum(log(pvalues_Pr_varied_primary_0.001))
pfinal_Pr_varied_primary_0.001<-pchisq(test.stat_Pr_varied_primary_0.001, df=2*length(pvalues_Pr_varied_primary_0.001), lower.tail=F)
pvalues_Pr_varied_primary_0.01 <- c(pvalue_Pr_varied_primary_0.01_1,pvalue_Pr_varied_primary_0.01_2,pvalue_Pr_varied_primary_0.01_3,pvalue_Pr_varied_primary_0.01_4,pvalue_Pr_varied_primary_0.01_5,pvalue_Pr_varied_primary_0.01_6,pvalue_Pr_varied_primary_0.01_7,pvalue_Pr_varied_primary_0.01_8,pvalue_Pr_varied_primary_0.01_9)
test.stat_Pr_varied_primary_0.01 <- -2*sum(log(pvalues_Pr_varied_primary_0.01))
pfinal_Pr_varied_primary_0.01<-pchisq(test.stat_Pr_varied_primary_0.01, df=2*length(pvalues_Pr_varied_primary_0.01), lower.tail=F)
pvalues_Pr_varied_primary_0.1 <- c(pvalue_Pr_varied_primary_0.1_1,pvalue_Pr_varied_primary_0.1_2,pvalue_Pr_varied_primary_0.1_3,pvalue_Pr_varied_primary_0.1_4,pvalue_Pr_varied_primary_0.1_5,pvalue_Pr_varied_primary_0.1_6,pvalue_Pr_varied_primary_0.1_7,pvalue_Pr_varied_primary_0.1_8,pvalue_Pr_varied_primary_0.1_9)
test.stat_Pr_varied_primary_0.1 <- -2*sum(log(pvalues_Pr_varied_primary_0.1))
pfinal_Pr_varied_primary_0.1<-pchisq(test.stat_Pr_varied_primary_0.1, df=2*length(pvalues_Pr_varied_primary_0.1), lower.tail=F)
pvalues_Pr_varied_primary_1.0 <- c(pvalue_Pr_varied_primary_1.0_1,pvalue_Pr_varied_primary_1.0_2,pvalue_Pr_varied_primary_1.0_3,pvalue_Pr_varied_primary_1.0_4,pvalue_Pr_varied_primary_1.0_5,pvalue_Pr_varied_primary_1.0_6,pvalue_Pr_varied_primary_1.0_7,pvalue_Pr_varied_primary_1.0_8,pvalue_Pr_varied_primary_1.0_9)
test.stat_Pr_varied_primary_1.0 <- -2*sum(log(pvalues_Pr_varied_primary_1.0))
pfinal_Pr_varied_primary_1.0<-pchisq(test.stat_Pr_varied_primary_1.0, df=2*length(pvalues_Pr_varied_primary_1.0), lower.tail=F)
pvalues_Pr_varied_primary_2.0 <- c(pvalue_Pr_varied_primary_2.0_1,pvalue_Pr_varied_primary_2.0_2,pvalue_Pr_varied_primary_2.0_3,pvalue_Pr_varied_primary_2.0_4,pvalue_Pr_varied_primary_2.0_5,pvalue_Pr_varied_primary_2.0_6,pvalue_Pr_varied_primary_2.0_8,pvalue_Pr_varied_primary_2.0_9)
test.stat_Pr_varied_primary_2.0 <- -2*sum(log(pvalues_Pr_varied_primary_2.0))
pfinal_Pr_varied_primary_2.0<-pchisq(test.stat_Pr_varied_primary_2.0, df=2*length(pvalues_Pr_varied_primary_2.0), lower.tail=F)
pvalues_Pr_varied_primary_0.001_table <- c(pvalue_Pr_varied_primary_0.001_1, pvalue_Pr_varied_primary_0.001_2, pvalue_Pr_varied_primary_0.001_3, pvalue_Pr_varied_primary_0.001_4, pvalue_Pr_varied_primary_0.001_5, pvalue_Pr_varied_primary_0.001_6, pvalue_Pr_varied_primary_0.001_7, pvalue_Pr_varied_primary_0.001_8, pvalue_Pr_varied_primary_0.001_9)
p_Pr_varied_primary_0.001<-scales::pvalue(pvalues_Pr_varied_primary_0.001_table)
pvalues_Pr_varied_primary_0.01_table <- c(pvalue_Pr_varied_primary_0.01_1,pvalue_Pr_varied_primary_0.01_2,pvalue_Pr_varied_primary_0.01_3,pvalue_Pr_varied_primary_0.01_4,pvalue_Pr_varied_primary_0.01_5,pvalue_Pr_varied_primary_0.01_6,pvalue_Pr_varied_primary_0.01_7,pvalue_Pr_varied_primary_0.01_8,pvalue_Pr_varied_primary_0.01_9)
p_Pr_varied_primary_0.01<-scales::pvalue(pvalues_Pr_varied_primary_0.01_table)
pvalues_Pr_varied_primary_0.1_table <- c(pvalue_Pr_varied_primary_0.1_1,pvalue_Pr_varied_primary_0.1_2, pvalue_Pr_varied_primary_0.1_3, pvalue_Pr_varied_primary_0.1_4, pvalue_Pr_varied_primary_0.1_5, pvalue_Pr_varied_primary_0.1_6, pvalue_Pr_varied_primary_0.1_7, pvalue_Pr_varied_primary_0.1_8, pvalue_Pr_varied_primary_0.1_9)
p_Pr_varied_primary_0.1<-scales::pvalue(pvalues_Pr_varied_primary_0.1_table)
pvalues_Pr_varied_primary_1.0_table <- c(pvalue_Pr_varied_primary_1.0_1,pvalue_Pr_varied_primary_1.0_2, pvalue_Pr_varied_primary_1.0_3, pvalue_Pr_varied_primary_1.0_4, pvalue_Pr_varied_primary_1.0_5, pvalue_Pr_varied_primary_1.0_6, pvalue_Pr_varied_primary_1.0_7, pvalue_Pr_varied_primary_1.0_8, pvalue_Pr_varied_primary_1.0_9)
p_Pr_varied_primary_1.0<-scales::pvalue(pvalues_Pr_varied_primary_1.0_table)
pvalues_Pr_varied_primary_2.0_table <- c(pvalue_Pr_varied_primary_2.0_1,pvalue_Pr_varied_primary_2.0_2, pvalue_Pr_varied_primary_2.0_3, pvalue_Pr_varied_primary_2.0_4, pvalue_Pr_varied_primary_2.0_5, pvalue_Pr_varied_primary_2.0_6, NA, pvalue_Pr_varied_primary_2.0_8, pvalue_Pr_varied_primary_2.0_9)
p_Pr_varied_primary_2.0<-scales::pvalue(pvalues_Pr_varied_primary_2.0_table)
p_Pr_varied_primary_2.0 <- p_Pr_varied_primary_2.0 %>% replace_na("--")
pvalue_table_Pr_varied_primary <- data.frame("0.001" = p_Pr_varied_primary_0.001, "0.01" = p_Pr_varied_primary_0.01, "0.1" = p_Pr_varied_primary_0.1, "1.0" = p_Pr_varied_primary_1.0, "2.0" = p_Pr_varied_primary_2.0)
colnames(pvalue_table_Pr_varied_primary) <- c(0.001,0.01, 0.1, 1.0,2.0)
pfisher_Pr_varied_primary<- c(pfinal_Pr_varied_primary_0.001, pfinal_Pr_varied_primary_0.01, pfinal_Pr_varied_primary_0.1, pfinal_Pr_varied_primary_1.0, NA)
pfisher_Pr_varied_primary<-scales::pvalue(pfisher_Pr_varied_primary)
pfisher_Pr_varied_primary<- pfisher_Pr_varied_primary%>% replace_na("N.C.")
pvalue_table_Pr_varied_primary[nrow(pvalue_table_Pr_varied_primary) + 1,]<-pfisher_Pr_varied_primary
row.names(pvalue_table_Pr_varied_primary) <- c("Block 1 - October 18th, 2020", "Block 2 - October 20th, 2020", "Block 3 - October 29th, 2020","Block 4 - November 2nd, 2020", "Block 5 - November 7th, 2020", "Block 6 - November 28th, 2020","Block 7 - December 4th, 2020", "Block 8 - June 17th, 2021", "Block 9 - June 18th, 2021", " Final P-value ")
pvalue_table_Pr_varied_primary %>%
kbl(caption = "<center><strong>P-values for log-rank tests on individual doses within each experiment<center><strong>") %>%
kable_classic(full_width = F, html_font = "Cambria", position = "center")%>%
add_header_above(c(" "=1, "Primary Infectious Dose (Abs600nm)" = 5)) %>%
pack_rows("Individual Blocks", 1, 9) %>%
pack_rows("Fisher's Combined", 10, 10)
|
Primary Infectious Dose (Abs600nm)
|
|||||
|---|---|---|---|---|---|
| 0.001 | 0.01 | 0.1 | 1 | 2 | |
| Individual Blocks | |||||
| Block 1 - October 18th, 2020 | 0.072 | 0.293 | <0.001 | 0.070 | 0.041 |
| Block 2 - October 20th, 2020 | 0.526 | 0.211 | <0.001 | 0.002 | 0.072 |
| Block 3 - October 29th, 2020 | 0.628 | 0.018 | <0.001 | <0.001 | 0.012 |
| Block 4 - November 2nd, 2020 | 0.942 | 0.009 | <0.001 | <0.001 | 0.407 |
| Block 5 - November 7th, 2020 | 0.693 | 0.332 | <0.001 | <0.001 | 0.010 |
| Block 6 - November 28th, 2020 | 0.064 | 0.913 | <0.001 | <0.001 | 0.033 |
| Block 7 - December 4th, 2020 | 0.968 | 0.007 | <0.001 | <0.001 | – |
| Block 8 - June 17th, 2021 | 0.726 | 0.063 | <0.001 | <0.001 | <0.001 |
| Block 9 - June 18th, 2021 | 0.238 | 0.001 | <0.001 | <0.001 | <0.001 |
| Fisher’s Combined | |||||
| Final P-value | 0.733 | <0.001 | <0.001 | <0.001 | N.C. |
Ps_varied_primary <- subset(Ps_varied_primary, secondary == "P.snee")
table(Ps_varied_primary$date, Ps_varied_primary$primary_dose)
##
## 0 0.001 0.01 0.1 1 2
## 3/10/2021 38 49 48 49 30 6
## 3/17/2021 37 50 46 42 22 0
## 3/31/2021 39 46 47 42 43 52
## 4/21/2021 35 44 46 45 30 49
## 4/28/2021 37 45 41 43 38 34
## 6/20/2021 27 48 49 48 47 60
table(Ps_varied_primary$primary_dose)
##
## 0 0.001 0.01 0.1 1 2
## 213 282 277 269 210 201
#png(file="Figure5B.png", width = 1200, height = 900, res = 130)
par(mar=c(6,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(Ps_varied_primary$hour, Ps_varied_primary$status) ~ Ps_varied_primary$condition), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=5, yaxt='n', xaxt='n', ylab="proportion alive", xlab="hours post-secondary infection", cex.lab=2.5, xlim=c(0,48), xaxs='i', ylim=c(0, 1), yaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
legend("bottomleft", bty="n" , c("secondary infection only", "chronic (0.001) & secondary infection", "chronic (0.01) & secondary infection", "chronic (0.1) & secondary infection", "chronic (1.0) & secondary infection", "chronic (2.0) & secondary infection"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1.5)
#dev.off()
Ps_varied_primary_block_1 <- subset(Ps_varied_primary, date == "3/10/2021")
table(Ps_varied_primary_block_1$primary_dose)
##
## 0 0.001 0.01 0.1 1 2
## 38 49 48 49 30 6
#subset by primary dose for pairwise comparison dose
Ps_varied_primary_0.001_1<-subset(Ps_varied_primary_block_1, primary_dose == "0.001"|primary_dose=="0")
table(Ps_varied_primary_0.001_1$primary_dose)
##
## 0 0.001
## 38 49
Ps_varied_primary_0.01_1<-subset(Ps_varied_primary_block_1, primary_dose=="0.01"|primary_dose=="0")
Ps_varied_primary_0.1_1<-subset(Ps_varied_primary_block_1, primary_dose=="0.1"|primary_dose=="0")
Ps_varied_primary_1.0_1<-subset(Ps_varied_primary_block_1, primary_dose=="1"|primary_dose=="0")
Ps_varied_primary_2.0_1<-subset(Ps_varied_primary_block_1, primary_dose=="2"|primary_dose=="0")
surv_Ps_varied_primary_0.001_1<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.001_1)
surv_Ps_varied_primary_0.001_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.001_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 38 38 30.3 1.98 4.76
## primary=S.m 49 47 54.7 1.10 4.76
##
## Chisq= 4.8 on 1 degrees of freedom, p= 0.03
surv_Ps_varied_primary_0.01_1<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.01_1)
surv_Ps_varied_primary_0.01_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.01_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 38 38 31.6 1.29 3.14
## primary=S.m 48 48 54.4 0.75 3.14
##
## Chisq= 3.1 on 1 degrees of freedom, p= 0.08
surv_Ps_varied_primary_0.1_1<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.1_1)
surv_Ps_varied_primary_0.1_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.1_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 38 38 22.1 11.40 22.7
## primary=S.m 49 49 64.9 3.89 22.7
##
## Chisq= 22.7 on 1 degrees of freedom, p= 2e-06
surv_Ps_varied_primary_1.0_1<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_1.0_1)
surv_Ps_varied_primary_1.0_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_1.0_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 38 38 36.7 0.0435 0.148
## primary=S.m 30 30 31.3 0.0511 0.148
##
## Chisq= 0.1 on 1 degrees of freedom, p= 0.7
surv_Ps_varied_primary_2.0_1<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_2.0_1)
surv_Ps_varied_primary_2.0_1
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_2.0_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 38 38 37.12 0.0208 0.234
## primary=S.m 6 5 5.88 0.1314 0.234
##
## Chisq= 0.2 on 1 degrees of freedom, p= 0.6
plot(survfit(Surv(Ps_varied_primary_block_1$hour, Ps_varied_primary_block_1$status) ~ Ps_varied_primary_block_1$condition), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Hours Post-Secondary Infection", cex.lab=1.3, xlim=c(10,48), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-P.snee 2.0", "S.m 0.001-P.snee 2.0", "S.m 0.01-P.snee 2.0", "S.m 0.1-P.snee 2.0", "S.m 1.0-P.snee 2.0", "S.m 2.0-P.snee 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)
Ps_varied_primary_block_2 <- subset(Ps_varied_primary, date == "3/17/2021")
table(Ps_varied_primary_block_2$primary_dose)
##
## 0 0.001 0.01 0.1 1
## 37 50 46 42 22
#subset by primary dose for pairwise comparison dose
Ps_varied_primary_0.001_2<-subset(Ps_varied_primary_block_2, primary_dose == "0.001"|primary_dose=="0")
table(Ps_varied_primary_0.001_2$primary_dose)
##
## 0 0.001
## 37 50
Ps_varied_primary_0.01_2<-subset(Ps_varied_primary_block_2, primary_dose=="0.01"|primary_dose=="0")
Ps_varied_primary_0.1_2<-subset(Ps_varied_primary_block_2, primary_dose=="0.1"|primary_dose=="0")
Ps_varied_primary_1.0_2<-subset(Ps_varied_primary_block_2, primary_dose=="1"|primary_dose=="0")
#Ps_varied_primary_2.0_2<-subset(Ps_varied_primary_block_2, primary_dose=="2"|primary_dose=="0")
surv_Ps_varied_primary_0.001_2<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.001_2)
surv_Ps_varied_primary_0.001_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.001_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 37 38.5 0.0568 0.16
## primary=S.m 50 49 47.5 0.0460 0.16
##
## Chisq= 0.2 on 1 degrees of freedom, p= 0.7
surv_Ps_varied_primary_0.01_2<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.01_2)
surv_Ps_varied_primary_0.01_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.01_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 37 30 1.638 3.92
## primary=S.m 46 44 51 0.963 3.92
##
## Chisq= 3.9 on 1 degrees of freedom, p= 0.05
surv_Ps_varied_primary_0.1_2<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.1_2)
surv_Ps_varied_primary_0.1_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.1_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 37 20.7 12.89 25.7
## primary=S.m 42 32 48.3 5.52 25.7
##
## Chisq= 25.7 on 1 degrees of freedom, p= 4e-07
surv_Ps_varied_primary_1.0_2<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_1.0_2)
surv_Ps_varied_primary_1.0_2
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_1.0_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 37 27.8 3.04 9.36
## primary=S.m 22 14 23.2 3.64 9.36
##
## Chisq= 9.4 on 1 degrees of freedom, p= 0.002
#surv_Ps_varied_primary_2.0_2<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_2.0_2)
#surv_Ps_varied_primary_2.0_2
plot(survfit(Surv(Ps_varied_primary_block_2$hour, Ps_varied_primary_block_2$status) ~ Ps_varied_primary_block_2$condition), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Hours Post-Secondary Infection", cex.lab=1.3, xlim=c(10,48), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-P.snee 2.0", "S.m 0.001-P.snee 2.0", "S.m 0.01-P.snee 2.0", "S.m 0.1-P.snee 2.0", "S.m 1.0-P.snee 2.0", "S.m 2.0-P.snee 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)
Ps_varied_primary_block_3 <- subset(Ps_varied_primary, date == "3/31/2021")
table(Ps_varied_primary_block_3$primary_dose)
##
## 0 0.001 0.01 0.1 1 2
## 39 46 47 42 43 52
#subset by primary dose for pairwise comparison dose
Ps_varied_primary_0.001_3<-subset(Ps_varied_primary_block_3, primary_dose == "0.001"|primary_dose=="0")
table(Ps_varied_primary_0.001_3$primary_dose)
##
## 0 0.001
## 39 46
Ps_varied_primary_0.01_3<-subset(Ps_varied_primary_block_3, primary_dose=="0.01"|primary_dose=="0")
Ps_varied_primary_0.1_3<-subset(Ps_varied_primary_block_3, primary_dose=="0.1"|primary_dose=="0")
Ps_varied_primary_1.0_3<-subset(Ps_varied_primary_block_3, primary_dose=="1"|primary_dose=="0")
Ps_varied_primary_2.0_3<-subset(Ps_varied_primary_block_3, primary_dose=="2"|primary_dose=="0")
surv_Ps_varied_primary_0.001_3<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.001_3)
surv_Ps_varied_primary_0.001_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.001_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 39 21 23 0.177 0.383
## primary=S.m 46 29 27 0.151 0.383
##
## Chisq= 0.4 on 1 degrees of freedom, p= 0.5
surv_Ps_varied_primary_0.01_3<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.01_3)
surv_Ps_varied_primary_0.01_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.01_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 39 21 20.4 0.0196 0.0399
## primary=S.m 47 26 26.6 0.0150 0.0399
##
## Chisq= 0 on 1 degrees of freedom, p= 0.8
surv_Ps_varied_primary_0.1_3<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.1_3)
surv_Ps_varied_primary_0.1_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.1_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 39 21 23.3 0.222 0.504
## primary=S.m 42 30 27.7 0.186 0.504
##
## Chisq= 0.5 on 1 degrees of freedom, p= 0.5
surv_Ps_varied_primary_1.0_3<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_1.0_3)
surv_Ps_varied_primary_1.0_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_1.0_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 39 21 12.3 6.08 11.3
## primary=S.m 43 8 16.7 4.50 11.3
##
## Chisq= 11.3 on 1 degrees of freedom, p= 8e-04
surv_Ps_varied_primary_2.0_3<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_2.0_3)
surv_Ps_varied_primary_2.0_3
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_2.0_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 39 21 10.7 9.96 17.2
## primary=S.m 52 7 17.3 6.15 17.2
##
## Chisq= 17.2 on 1 degrees of freedom, p= 3e-05
plot(survfit(Surv(Ps_varied_primary_block_3$hour, Ps_varied_primary_block_3$status) ~ Ps_varied_primary_block_3$condition), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Hours Post-Secondary Infection", cex.lab=1.3, xlim=c(10,48), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-P.snee 2.0", "S.m 0.001-P.snee 2.0", "S.m 0.01-P.snee 2.0", "S.m 0.1-P.snee 2.0", "S.m 1.0-P.snee 2.0", "S.m 2.0-P.snee 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)
Ps_varied_primary_block_4 <- subset(Ps_varied_primary, date == "4/21/2021")
table(Ps_varied_primary_block_4$primary_dose)
##
## 0 0.001 0.01 0.1 1 2
## 35 44 46 45 30 49
#subset by primary dose for pairwise comparison dose
Ps_varied_primary_0.001_4<-subset(Ps_varied_primary_block_4, primary_dose == "0.001"|primary_dose=="0")
table(Ps_varied_primary_0.001_4$primary_dose)
##
## 0 0.001
## 35 44
Ps_varied_primary_0.01_4<-subset(Ps_varied_primary_block_4, primary_dose=="0.01"|primary_dose=="0")
Ps_varied_primary_0.1_4<-subset(Ps_varied_primary_block_4, primary_dose=="0.1"|primary_dose=="0")
Ps_varied_primary_1.0_4<-subset(Ps_varied_primary_block_4, primary_dose=="1"|primary_dose=="0")
Ps_varied_primary_2.0_4<-subset(Ps_varied_primary_block_4, primary_dose=="2"|primary_dose=="0")
surv_Ps_varied_primary_0.001_4<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.001_4)
surv_Ps_varied_primary_0.001_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.001_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 35 35 26.4 2.83 6.39
## primary=S.m 44 44 52.6 1.42 6.39
##
## Chisq= 6.4 on 1 degrees of freedom, p= 0.01
surv_Ps_varied_primary_0.01_4<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.01_4)
surv_Ps_varied_primary_0.01_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.01_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 35 35 18.9 13.74 26
## primary=S.m 46 39 55.1 4.71 26
##
## Chisq= 26 on 1 degrees of freedom, p= 3e-07
surv_Ps_varied_primary_0.1_4<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.1_4)
surv_Ps_varied_primary_0.1_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.1_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 35 35 12.4 41.3 69.8
## primary=S.m 45 22 44.6 11.5 69.8
##
## Chisq= 69.8 on 1 degrees of freedom, p= <2e-16
surv_Ps_varied_primary_1.0_4<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_1.0_4)
surv_Ps_varied_primary_1.0_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_1.0_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 35 35 16 22.5 46.5
## primary=S.m 30 7 26 13.9 46.5
##
## Chisq= 46.5 on 1 degrees of freedom, p= 9e-12
surv_Ps_varied_primary_2.0_4<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_2.0_4)
surv_Ps_varied_primary_2.0_4
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_2.0_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 35 35 19.5 12.33 23
## primary=S.m 49 26 41.5 5.79 23
##
## Chisq= 23 on 1 degrees of freedom, p= 2e-06
plot(survfit(Surv(Ps_varied_primary_block_4$hour, Ps_varied_primary_block_4$status) ~ Ps_varied_primary_block_4$condition), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Hours Post-Secondary Infection", cex.lab=1.3, xlim=c(10,48), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-P.snee 2.0", "S.m 0.001-P.snee 2.0", "S.m 0.01-P.snee 2.0", "S.m 0.1-P.snee 2.0", "S.m 1.0-P.snee 2.0", "S.m 2.0-P.snee 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)
Ps_varied_primary_block_5 <- subset(Ps_varied_primary, date == "4/28/2021")
table(Ps_varied_primary_block_5$primary_dose)
##
## 0 0.001 0.01 0.1 1 2
## 37 45 41 43 38 34
#subset by primary dose for pairwise comparison dose
Ps_varied_primary_0.001_5<-subset(Ps_varied_primary_block_5, primary_dose == "0.001"|primary_dose=="0")
table(Ps_varied_primary_0.001_5$primary_dose)
##
## 0 0.001
## 37 45
Ps_varied_primary_0.01_5<-subset(Ps_varied_primary_block_5, primary_dose=="0.01"|primary_dose=="0")
Ps_varied_primary_0.1_5<-subset(Ps_varied_primary_block_5, primary_dose=="0.1"|primary_dose=="0")
Ps_varied_primary_1.0_5<-subset(Ps_varied_primary_block_5, primary_dose=="1"|primary_dose=="0")
Ps_varied_primary_2.0_5<-subset(Ps_varied_primary_block_5, primary_dose=="2"|primary_dose=="0")
surv_Ps_varied_primary_0.001_5<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.001_5)
surv_Ps_varied_primary_0.001_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.001_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 37 37.2 0.000842 0.00276
## primary=S.m 45 45 44.8 0.000699 0.00276
##
## Chisq= 0 on 1 degrees of freedom, p= 1
surv_Ps_varied_primary_0.01_5<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.01_5)
surv_Ps_varied_primary_0.01_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.01_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 37 36.4 0.00838 0.0276
## primary=S.m 41 41 41.6 0.00735 0.0276
##
## Chisq= 0 on 1 degrees of freedom, p= 0.9
surv_Ps_varied_primary_0.1_5<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.1_5)
surv_Ps_varied_primary_0.1_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.1_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 37 28.3 2.68 6.9
## primary=S.m 43 43 51.7 1.47 6.9
##
## Chisq= 6.9 on 1 degrees of freedom, p= 0.009
surv_Ps_varied_primary_1.0_5<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_1.0_5)
surv_Ps_varied_primary_1.0_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_1.0_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 37 17.9 20.24 40.2
## primary=S.m 38 34 53.1 6.85 40.2
##
## Chisq= 40.2 on 1 degrees of freedom, p= 2e-10
surv_Ps_varied_primary_2.0_5<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_2.0_5)
surv_Ps_varied_primary_2.0_5
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_2.0_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 37 18.8 17.71 37.4
## primary=S.m 34 30 48.2 6.89 37.4
##
## Chisq= 37.4 on 1 degrees of freedom, p= 9e-10
plot(survfit(Surv(Ps_varied_primary_block_5$hour, Ps_varied_primary_block_5$status) ~ Ps_varied_primary_block_5$condition), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Hours Post-Secondary Infection", cex.lab=1.3, xlim=c(10,48), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-P.snee 2.0", "S.m 0.001-P.snee 2.0", "S.m 0.01-P.snee 2.0", "S.m 0.1-P.snee 2.0", "S.m 1.0-P.snee 2.0", "S.m 2.0-P.snee 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)
Ps_varied_primary_block_6 <- subset(Ps_varied_primary, date == "6/20/2021")
table(Ps_varied_primary_block_6$primary_dose)
##
## 0 0.001 0.01 0.1 1 2
## 27 48 49 48 47 60
#subset by primary dose for pairwise comparison dose
Ps_varied_primary_0.001_6<-subset(Ps_varied_primary_block_6, primary_dose == "0.001"|primary_dose=="0")
table(Ps_varied_primary_0.001_6$primary_dose)
##
## 0 0.001
## 27 48
Ps_varied_primary_0.01_6<-subset(Ps_varied_primary_block_6, primary_dose=="0.01"|primary_dose=="0")
Ps_varied_primary_0.1_6<-subset(Ps_varied_primary_block_6, primary_dose=="0.1"|primary_dose=="0")
Ps_varied_primary_1.0_6<-subset(Ps_varied_primary_block_6, primary_dose=="1"|primary_dose=="0")
Ps_varied_primary_2.0_6<-subset(Ps_varied_primary_block_6, primary_dose=="2"|primary_dose=="0")
surv_Ps_varied_primary_0.001_6<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.001_6)
surv_Ps_varied_primary_0.001_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.001_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 27 27 24.9 0.1730 0.436
## primary=S.m 48 48 50.1 0.0861 0.436
##
## Chisq= 0.4 on 1 degrees of freedom, p= 0.5
surv_Ps_varied_primary_0.01_6<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.01_6)
surv_Ps_varied_primary_0.01_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.01_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 27 27 20.2 2.262 4.59
## primary=S.m 49 46 52.8 0.868 4.59
##
## Chisq= 4.6 on 1 degrees of freedom, p= 0.03
surv_Ps_varied_primary_0.1_6<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.1_6)
surv_Ps_varied_primary_0.1_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.1_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 27 27 10.5 25.69 40.5
## primary=S.m 48 31 47.5 5.71 40.5
##
## Chisq= 40.5 on 1 degrees of freedom, p= 2e-10
surv_Ps_varied_primary_1.0_6<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_1.0_6)
surv_Ps_varied_primary_1.0_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_1.0_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 27 27 7.82 47.01 71.3
## primary=S.m 47 19 38.18 9.63 71.3
##
## Chisq= 71.3 on 1 degrees of freedom, p= <2e-16
surv_Ps_varied_primary_2.0_6<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_2.0_6)
surv_Ps_varied_primary_2.0_6
## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_2.0_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 27 27 14.1 11.88 19
## primary=S.m 60 37 49.9 3.35 19
##
## Chisq= 19 on 1 degrees of freedom, p= 1e-05
plot(survfit(Surv(Ps_varied_primary_block_6$hour, Ps_varied_primary_block_6$status) ~ Ps_varied_primary_block_6$condition), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Hours Post-Secondary Infection", cex.lab=1.3, xlim=c(10,48), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-P.snee 2.0", "S.m 0.001-P.snee 2.0", "S.m 0.01-P.snee 2.0", "S.m 0.1-P.snee 2.0", "S.m 1.0-P.snee 2.0", "S.m 2.0-P.snee 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)
# Grab p-values from each individual logrank test (results shown above in each block)
pvalue_Ps_varied_primary_0.001_1<-pchisq(surv_Ps_varied_primary_0.001_1$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.01_1<-pchisq(surv_Ps_varied_primary_0.01_1$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.1_1<-pchisq(surv_Ps_varied_primary_0.1_1$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_1.0_1<-pchisq(surv_Ps_varied_primary_1.0_1$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_2.0_1<-pchisq(surv_Ps_varied_primary_2.0_1$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.001_2<-pchisq(surv_Ps_varied_primary_0.001_2$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.01_2<-pchisq(surv_Ps_varied_primary_0.01_2$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.1_2<-pchisq(surv_Ps_varied_primary_0.1_2$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_1.0_2<-pchisq(surv_Ps_varied_primary_1.0_2$chisq, df=1, lower.tail=F)
#pvalue_Ps_varied_primary_2.0_2<-pchisq(surv_Ps_varied_primary_2.0_2$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.001_3<-pchisq(surv_Ps_varied_primary_0.001_3$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.01_3<-pchisq(surv_Ps_varied_primary_0.01_3$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.1_3<-pchisq(surv_Ps_varied_primary_0.1_3$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_1.0_3<-pchisq(surv_Ps_varied_primary_1.0_3$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_2.0_3<-pchisq(surv_Ps_varied_primary_2.0_3$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.001_4<-pchisq(surv_Ps_varied_primary_0.001_4$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.01_4<-pchisq(surv_Ps_varied_primary_0.01_4$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.1_4<-pchisq(surv_Ps_varied_primary_0.1_4$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_1.0_4<-pchisq(surv_Ps_varied_primary_1.0_4$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_2.0_4<-pchisq(surv_Ps_varied_primary_2.0_4$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.001_5<-pchisq(surv_Ps_varied_primary_0.001_5$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.01_5<-pchisq(surv_Ps_varied_primary_0.01_5$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.1_5<-pchisq(surv_Ps_varied_primary_0.1_5$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_1.0_5<-pchisq(surv_Ps_varied_primary_1.0_5$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_2.0_5<-pchisq(surv_Ps_varied_primary_2.0_5$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.001_6<-pchisq(surv_Ps_varied_primary_0.001_6$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.01_6<-pchisq(surv_Ps_varied_primary_0.01_6$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.1_6<-pchisq(surv_Ps_varied_primary_0.1_6$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_1.0_6<-pchisq(surv_Ps_varied_primary_1.0_6$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_2.0_6<-pchisq(surv_Ps_varied_primary_2.0_6$chisq, df=1, lower.tail=F)
# Use all p-values from a given condition to get Fisher's combined test p-value
pvalues_Ps_varied_primary_0.001 <- c(pvalue_Ps_varied_primary_0.001_2 ,pvalue_Ps_varied_primary_0.001_3, pvalue_Ps_varied_primary_0.001_4, pvalue_Ps_varied_primary_0.001_5, pvalue_Ps_varied_primary_0.001_6)
test.stat_Ps_varied_primary_0.001 <- -2*sum(log(pvalues_Ps_varied_primary_0.001))
pfinal_Ps_varied_primary_0.001<-pchisq(test.stat_Ps_varied_primary_0.001, df=2*length(pvalues_Ps_varied_primary_0.001), lower.tail=F)
pvalues_Ps_varied_primary_0.01 <- c(pvalue_Ps_varied_primary_0.01_1,pvalue_Ps_varied_primary_0.01_2,pvalue_Ps_varied_primary_0.01_3,pvalue_Ps_varied_primary_0.01_4,pvalue_Ps_varied_primary_0.01_5,pvalue_Ps_varied_primary_0.01_6)
test.stat_Ps_varied_primary_0.01 <- -2*sum(log(pvalues_Ps_varied_primary_0.01))
pfinal_Ps_varied_primary_0.01<-pchisq(test.stat_Ps_varied_primary_0.01, df=2*length(pvalues_Ps_varied_primary_0.01), lower.tail=F)
pvalues_Ps_varied_primary_0.1 <- c(pvalue_Ps_varied_primary_0.1_1,pvalue_Ps_varied_primary_0.1_2,pvalue_Ps_varied_primary_0.1_3,pvalue_Ps_varied_primary_0.1_4,pvalue_Ps_varied_primary_0.1_5,pvalue_Ps_varied_primary_0.1_6)
test.stat_Ps_varied_primary_0.1 <- -2*sum(log(pvalues_Ps_varied_primary_0.1))
pfinal_Ps_varied_primary_0.1<-pchisq(test.stat_Ps_varied_primary_0.1, df=2*length(pvalues_Ps_varied_primary_0.1), lower.tail=F)
pvalues_Ps_varied_primary_1.0 <- c(pvalue_Ps_varied_primary_1.0_1,pvalue_Ps_varied_primary_1.0_2,pvalue_Ps_varied_primary_1.0_3,pvalue_Ps_varied_primary_1.0_4,pvalue_Ps_varied_primary_1.0_5,pvalue_Ps_varied_primary_1.0_6)
test.stat_Ps_varied_primary_1.0 <- -2*sum(log(pvalues_Ps_varied_primary_1.0))
pfinal_Ps_varied_primary_1.0<-pchisq(test.stat_Ps_varied_primary_1.0, df=2*length(pvalues_Ps_varied_primary_1.0), lower.tail=F)
pvalues_Ps_varied_primary_2.0 <- c(pvalue_Ps_varied_primary_2.0_1,pvalue_Ps_varied_primary_2.0_3,pvalue_Ps_varied_primary_2.0_4,pvalue_Ps_varied_primary_2.0_5,pvalue_Ps_varied_primary_2.0_6)
test.stat_Ps_varied_primary_2.0 <- -2*sum(log(pvalues_Ps_varied_primary_2.0))
pfinal_Ps_varied_primary_2.0<-pchisq(test.stat_Ps_varied_primary_2.0, df=2*length(pvalues_Ps_varied_primary_2.0), lower.tail=F)
pvalues_Ps_varied_primary_0.001_table <- c(pvalue_Ps_varied_primary_0.001_1, pvalue_Ps_varied_primary_0.001_2, pvalue_Ps_varied_primary_0.001_3, pvalue_Ps_varied_primary_0.001_4, pvalue_Ps_varied_primary_0.001_5, pvalue_Ps_varied_primary_0.001_6)
p_Ps_varied_primary_0.001<-scales::pvalue(pvalues_Ps_varied_primary_0.001_table)
pvalues_Ps_varied_primary_0.01_table <- c(pvalue_Ps_varied_primary_0.01_1,pvalue_Ps_varied_primary_0.01_2,pvalue_Ps_varied_primary_0.01_3,pvalue_Ps_varied_primary_0.01_4,pvalue_Ps_varied_primary_0.01_5,pvalue_Ps_varied_primary_0.01_6)
p_Ps_varied_primary_0.01<-scales::pvalue(pvalues_Ps_varied_primary_0.01_table)
pvalues_Ps_varied_primary_0.1_table <- c(pvalue_Ps_varied_primary_0.1_1,pvalue_Ps_varied_primary_0.1_2, pvalue_Ps_varied_primary_0.1_3, pvalue_Ps_varied_primary_0.1_4, pvalue_Ps_varied_primary_0.1_5, pvalue_Ps_varied_primary_0.1_6)
p_Ps_varied_primary_0.1<-scales::pvalue(pvalues_Ps_varied_primary_0.1_table)
pvalues_Ps_varied_primary_1.0_table <- c(pvalue_Ps_varied_primary_1.0_1,pvalue_Ps_varied_primary_1.0_2, pvalue_Ps_varied_primary_1.0_3, pvalue_Ps_varied_primary_1.0_4, pvalue_Ps_varied_primary_1.0_5, pvalue_Ps_varied_primary_1.0_6)
p_Ps_varied_primary_1.0<-scales::pvalue(pvalues_Ps_varied_primary_1.0_table)
pvalues_Ps_varied_primary_2.0_table <- c(pvalue_Ps_varied_primary_2.0_1,NA, pvalue_Ps_varied_primary_2.0_3, pvalue_Ps_varied_primary_2.0_4, pvalue_Ps_varied_primary_2.0_5, pvalue_Ps_varied_primary_2.0_6)
p_Ps_varied_primary_2.0<-scales::pvalue(pvalues_Ps_varied_primary_2.0_table)
p_Ps_varied_primary_2.0 <- p_Ps_varied_primary_2.0 %>% replace_na("--")
pvalue_table_Ps_varied_primary <- data.frame("0.001" = p_Ps_varied_primary_0.001, "0.01" = p_Ps_varied_primary_0.01, "0.1" = p_Ps_varied_primary_0.1, "1.0" = p_Ps_varied_primary_1.0, "2.0" = p_Ps_varied_primary_2.0)
colnames(pvalue_table_Ps_varied_primary) <- c(0.001,0.01, 0.1, 1.0,2.0)
pfisher_Ps_varied_primary<- c(pfinal_Ps_varied_primary_0.001, pfinal_Ps_varied_primary_0.01, pfinal_Ps_varied_primary_0.1, pfinal_Ps_varied_primary_1.0, pfinal_Ps_varied_primary_2.0)
pfisher_Ps_varied_primary<-scales::pvalue(pfisher_Ps_varied_primary)
pvalue_table_Ps_varied_primary[nrow(pvalue_table_Ps_varied_primary) + 1,]<-pfisher_Ps_varied_primary
row.names(pvalue_table_Ps_varied_primary) <- c("Block 1 - March 10th, 2021", "Block 2 - March 17th, 2021", "Block 3 - March 31st, 2021","Block 4 - April 21t, 2021", "Block 5 - April 28th, 2021", "Block 6 - June 20th, 2021", " Final P-value ")
pvalue_table_Ps_varied_primary %>%
kbl(caption = "<center><strong>P-values for log-rank tests on individual doses within each experiment<center><strong>") %>%
kable_classic(full_width = F, html_font = "Cambria", position = "center")%>%
add_header_above(c(" "=1, "Primary Infectious Dose (Abs600nm)" = 5)) %>%
pack_rows("Individual Blocks", 1, 6) %>%
pack_rows("Fisher's Combined", 7, 7)
|
Primary Infectious Dose (Abs600nm)
|
|||||
|---|---|---|---|---|---|
| 0.001 | 0.01 | 0.1 | 1 | 2 | |
| Individual Blocks | |||||
| Block 1 - March 10th, 2021 | 0.029 | 0.076 | <0.001 | 0.700 | 0.628 |
| Block 2 - March 17th, 2021 | 0.689 | 0.048 | <0.001 | 0.002 | – |
| Block 3 - March 31st, 2021 | 0.536 | 0.842 | 0.478 | <0.001 | <0.001 |
| Block 4 - April 21t, 2021 | 0.011 | <0.001 | <0.001 | <0.001 | <0.001 |
| Block 5 - April 28th, 2021 | 0.958 | 0.868 | 0.009 | <0.001 | <0.001 |
| Block 6 - June 20th, 2021 | 0.509 | 0.032 | <0.001 | <0.001 | <0.001 |
| Fisher’s Combined | |||||
| Final P-value | 0.261 | <0.001 | <0.001 | <0.001 | <0.001 |
q <- ggdraw() +
draw_image("Figure5A.png")
r <- ggdraw() +
draw_image("Figure5B.png")
multiplot5 <- q/r
v<-multiplot5
#plot_annotation(tag_levels = 'A')
ggsave("Fig.5.Primary.pdf", v, device = "pdf")
v